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EXECUTIVE  SUMMARY 


OBJECTIVE 

Develop  a  computer  code  to  predict  sea  radiance  (brightness). 

APPROACH 

Sea  radiance  is  modeled  by  combining  the  methods  of  geometrical  optics  with  the  Cox-Munk  sta¬ 
tistical  description  of  ocean  capillary  waves.  The  model  is  incorporated  into  the  atmospheric  transmit¬ 
tance/radiance  code  MODTRAN2  to  provide  numerical  sea  radiance  predictions. 

In  this  model  each  individual  capillary  wave  facet  is  allowed  to  reflect  the  sky  or  sun  and  emit  ther¬ 
mal  radiation.  The  total  radiance  from  the  sea  is  obtained  by  applying  the  proper  statistical  weight  to 
each  facet  and  integrating  over  all  facets  within  the  observer’s  field-of-view. 

RESULTS 

The  modified  MODTRAN2  code,  called  SeaRad,  calculates  sea  radiance  for  any  viewing  geometry 
in  a  spectral  range  from  52.63  cm'*  to  25000  cm'*.  Typical  execution  speeds  are  approximately  10  s 
per  pixel  on  a  Pentium/90  MHz  personal  computer.  Preliminary  comparisons  show  that  SeaRad 
agrees  to  within  several  degrees  Celsius  (°C)  with  actual  sea  radiance  measurements  in  the  mid-wave 
and  long-wave  infrared  bands. 
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1.  INTRODUCTION 


SeaRad  is  a  FORTRAN  computer  code  that  predicts  the  radiance  (brightness)  of  the  ocean  surface. 
SeaRad  is  based  on  the  Cox-Munk  statistical  model  (Cox  and  Munk,  1954,  1956)  for  wind-driven 
capillary  wave  facets.  An  individual  facet  is  chosen  and  assigned  a  specific  slope  with  respect  to  the 
local  horizon.  The  facet  is  allowed  to  reflect  the  sky  and  sun  and  emit  thermal  black  body  radiation 
toward  an  observer.  The  total  radiance  is  obtained  by  applying  the  proper  statistical  weight  to  the 
facet  and  integrating  over  all  facets  within  the  observer’s  field-of-view. 

SeaRad  is  valid  for  a  spectral  range  extending  from  the  visible  to  the  far  infrared.  Preliminary 
comparisons  show  that  SeaRad  agrees  to  within  several  °C  with  actual  sea  radiance  measurements  in 
the  mid- wave  and  long-wave  infrared  bands. 

In  its  current  form,  SeaRad  is  a  self-contained,  DOS-compatible  program  that  runs  on  a  personal 
computer  and  computes  radiance  for  a  single  pixel  (rather  than  an  entire  image).  It  is  a  modified  ver¬ 
sion  of  the  Air  Force  program  MODTRAN2  (Berk,  et  al.,  1989;  Kneizys,  et  al.,  1988)  that  computes 
atmospheric  transmittance  and  radiance.  SeaRad  operates  exactly  like  the  original  MODTRAN2 
code!  except  that  a  new  logical  parameter,  “SeaS witch”,  is  required  in  the  input  file.  Sun  glint  is 
included  in  the  sea  radiance  prediction  provided  that  the  user  has  chosen  to  execute  SeaRad  in 
radiance  mode  with  solar  scattered  radiance  included  (lEMSCT  =  2). 

2.  HARDWARE  CONSIDERATIONS 

The  size  of  the  FORTRAN  source  code  is  1.8  MB.  When  assembled  by  version  5.01  of  the  Lahey 
F77L/EM-32  DOS  compiler,  the  size  of  the  executable  code  is  0.8  MB.  When  run^  on  a  Pentium/90 
at  low  spectral  resolution  (LOWTRAN7)  in  multiple  scattering  mode,  execution  times  are  4  s  for  a 
typical  thermal  long-wave  case  (830  to  1250  cm'^  in  21  spectral  steps)  and  17  s  for  a  typical  solar 
mid-wave  case  (2000  to  3340  cm'^  in  67  spectral  steps).  Source  and  executable  codes  are  available 
on  disk  through  correspondence  with  the  author. 

3.  AN  EXAMPLE 

This  section  provides  an  example  of  how  SeaRad  is  used  to  predict  radiance  of  the  ocean  surface. 
An  input  file  called  “Tape5Rad.Std,”  shown  on  page  A-2,  employs  a  1976  U.  S.  standard  atmosphere 
to  calculate  ocean  radiance  observed  at  a  zenith  angle  of  100°  (a  depression  angle  of  10°)  from  a 
height  of  23  m.  The  Navy  aerosol  model  is  used.  The  calculation  is  done  at  low  spectral  resolution 
(LOWTRAN7)  for  a  single  wave  number  (945  cm’^)  in  the  long-wave  band. 

With  this  file  present,  the  following  three  DOS  commands  will  calculate  ocean  radiance  and  dis¬ 
play  results: 

Copy  Tape5Rad.Std  Tape5 

SeaRad 

Type  Out 


1.  This  report  assumes  that  the  reader  is  familiar  with  MODTRAN2  operation. 

2.  The  compiler  requires  the  Lahey/Phar  Lap  386  DOS  Extender  program  (0.2  MB)  to  run  on  a  personal  computer. 
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These  commands^  produce  the  output  file  “Out”  (Appendix  A,  page  3).  Band-integrated  radiance 
values  in  W  m’^  sr^  are  listed  at  the  end  of  the  output  file  for  each  of  four  contributions  to  ocean 
radiance:  path  to  footprint,  sea  emission,  sky  reflection,  and  sun  glint.  (In  fact,  no  sun  glint  has  been 
calculated  in  this  instance  since  the  input  file  specifies  lEMSCT  =  1  rather  than  lEMSCT  =  2.)  Please 
note  that  the  parameter  “TBOUND”  in  the  input  file  has  been  reinterpreted  by  SeaRad  as  the  sea 
temperature. 

The  input  file  shown  in  Appendix  A  on  page  2  contains  two  new  parameters  at  the  end  of  the  third 
line:  “90.000”  and  “T”.  These  will  be  discussed  in  reverse  order  of  their  appearance. 

The  “T”,  which  may  appear  anywhere  in  columns  76  through  80  of  the  third  line  of  the  input  file 
(at  the  end’of  Card  3),  is  a  new  logical  parameter  “SeaSwitch”.  It  is  required;  that  is,  a  fatal  error  will 
be  generated  if  it  is  not  present  in  the  input  file.  “SeaSwitch”  controls  the  sea  radiance  calculation. 
When  “SeaSwitch”  is  equal  to  “T”,  the  sea  radiance  calculation  will  be  allowed  provided  certain 
other  conditions  are  met.  When  “SeaSwitch”  is  equal  to  “F’,  the  sea  radiance  calculation  will  be  pre¬ 
vented  under  all  conditions  and  the  program  will  execute  as  originally  released  by  the  Air  Force. 

The  “90.000”,  which  may  appear  anywhere  in  columns  66  through  75  of  the  third  line  of  the  input 
file  (near  the  end  of  Card  3),  is  a  new  floating  point  parameter,  “Psi”.  It  is  optional;  that  is,  the  pro¬ 
gram  will  run  whether  this  parameter  is  included  in  the  input  file  or  not.  Psi  is  the  azimuth  of  the 
upwind  direction'*  measured  from  the  line-of-sight  in  degrees  positive  East  of  North.  If  it  is  omitted 
(if  the  field  is  blank),  and  if  all  conditions  for  a  sea  radiance  calculation  are  met,  that  calculation  will 
proceed  under  the  assumption  that  the  value  of  “Psi”  is  zero,  meaning  that  the  observer  is  looking 
directly  into  the  wind.  For  the  input  file  in  Appendix  A,  “Psi”  is  90°,  meaning  that  the  wind  is  blow¬ 
ing  from  right  to  left,  perpendicular  to  the  direction  of  observation. 

The  modified  version  of  Card  3  used  by  SeaRad  is: 

HI,  H2,  ANGLE,  RANGE,  BETA,  RO,  LEN,  Psi,  SeaSwitch 
Format  (6F10.3, 15,  F10.3,  L5) 

4.  THE  MODEL 

The  primary  assumption  of  the  model  is  that  the  strength  of  interaction  between  an  optical  ray  and 
a  capillary  wave  facet  is  given  by  the  facet  area  projected  normal  to  the  ray.  A  feature  (Zeisse,  1994. 
1995)  of  the  equations  contained  in  SeaRad  is  that  they  predict  a  finite  horizon  radiance.  SeaRad 
does  not  include  multiple  reflections,  shadowing,  or  gravity  waves.  Polarization  is  ignored. 

The  model  computes  four  contributions  to  sea  radiance.  Each  of  the  four  contributions  is  shown  in 
figure  1.  (For  purposes  of  clarity,  only  two  dimensions  have  been  used  in  figure  1;  however,  all  three 
dimensions  are  used  in  the  actual  calculation.) 

The  first  contribution  is  path  radiance,  shown  at  the  top  of  figure  1 .  The  footprint  of  a  single  pixel 
in  an  image  of  the  sea  is  indicated  by  the  wavy  line.  The  footprint  is  observed  by  a  receiver  at  the  end 


3.  The  time  for  this  particular  test  case  was  3  s  on  a  486/50  MHz  personal  computer. 

4.  This  information  is  required  because  the  Cox-Munk  capillary  wave  slope  statistics  are  different  in  the  upwind  and 
crosswind  directions. 
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of  a  ray  whose  zenith  angle  at  the  footprint  is  Qy.  Let  l^path  designate  the  spectral  radiance  in  W 
sr'  (cnT^)'i  along  the  path^  from  the  footprint  to  the  receiver. 

The  second  contribution  is  reflected  sky  radiance.  Spectral  radiance  from  a  portion  of  the  sky 
arrives  at  the  footprint  along  a  ray  whose  zenith  angle  there  is  d^-  The  footprint  contains  wave  facets 
of  different  slopes,  many  that  reflect  the  incoming  sky  radiance  away  from  the  receiver  toward  other 
parts  of  the  sky.  These  facets  are  ignored.  However,  the  footprint  will  contain  some  facets  whose 
slope  is  correct  for  reflecting  the  incoming  sky  radiance  toward  the  receiver  along  the  path  defined 
by  the  zenith  angle  Qy.  These  facets  are  retained.  The  contributions  from  all  portions  of  the  sky  are 
summed  together  after  specular  reflection  by  the  appropriate  facets  within  the  footprint,  and  the  sum 
leaving  the  footprint  at  zenith  angle  Qy  is  designated  During  its  path  to  the  receiver,  the  reflected 
sky  radiance  is  attenuated  by  the  path  transmission  Xpath- 

The  third  contribution  is  reflected  solar  radiance,  sun  glint.  The  calculation  is  analogous  to  the  cal¬ 
culation  of  sky  radiance.  Spectral  radiance  Nq  from  the  solar  center  arrives  at  the  footprint  along  a 
path  whose  zenith  angle  there  is  Qp-  Within  the  footprint,  most  facets  deflect  the  solar  ray  away  from 
the  receiver  and  are  rejected,  but  some  facets  are  retained  because  they  deflect  the  ray  specularly 
toward  the  receiver  along  a  path  with  zenith  angle  Qy.  N^un  is  the  spectral  radiance  leaving  the  foot¬ 
print  after  summation  over  rays  arriving  from  all  portions  of  the  solar  disk.  The  reflected  solar 
radiance  is  also  attenuated  by  the  path  transmission  Xpath  before  final  reception. 

The  fourth  contribution  is  thermal  black  body  emission.  Each  facet  emits  a  spectral  radiance  N^t, 
given  by  Planck’s  equation  for  a  black  body  whose  temperature  is  equal  to  the  value  of  “TBOUND” 
in  the  input  file.  The  spectral  emissivity  of  a  given  facet  in  the  direction  of  the  receiver  is  specified 
by  the  slope  of  that  facet  and  the  value  of  dy.  Nsea  is  the  thermal  spectral  radiance  leaving  the  foot¬ 
print  for  the  receiver  after  summation  over  all  facets  within  the  footprint.  As  before,  N^ea  is  atte¬ 
nuated  by  path  transmission  after  leaving  the  footprint. 

Throughout  figure  1,  the  symbol  p  represents  the  spectral  reflectivity  of  sea  water,  which  is 
required  for  the  second  and  third  contributions  since  they  are  governed  by  the  process  of  optical 
reflection.  On  the  other  hand,  the  fourth  contribution  is  governed  by  the  process  of  optical  emission. 
Fortunately,  by  application  of  Kirchoff’s  Law  to  an  opaque  medium,  sea  water,  the  emissivity  is 
given  by  one  minus  the  reflectivity.  The  reflectivity  is  calculated  from  Fresnel’s  equations  (Stratton, 
1941)  with  a  complex  optical  index  taken  from  the  literature  (Hale  &  Querry,  1973;  Querry,  et  al., 
1977).  These  data  for  the  index,  available  between  52.63  cm’’  and  25000  cm-^  set  the  spectral  range 
of  SeaRad. 


The  total  spectral  radiance  N  received  at  wave  number  v  (cm'*)  is  given  by 


Mv)  =  Va(v)/(v)  + 


Nskyiv) 


(1) 


where /(v)  has  been  introduced  to  represent  the  spectral  responsivity  of  the  receiver. 

The  design  of  SeaRad  is  such  that  path  {Npath->  T^path’)  ^nd  source  (Ns,  Np,  N^h)  values  are  taken 
from  the  original  MODTRAN2  while  Fresnel  reflection  (p)  and  slope  integrated  values  (NsPy,  Nsm^ 
Nsea)  are  introduced  in  new  subroutines.  Integration  of  (1)  over  the  wave  number  band  specified  in 
the  input  file  is  carried  out  in  a  modification  of  subroutine  “TRANS”  to  produce  the  band-integrated 
values  for  sea  radiance  given  in  the  output  file. 

5.  In  this  report,  the  word  path  refers  to  only  the  optical  path  between  the  footprint  and  the  receiver. 
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5.  THE  COORDINATE  SYSTEM 


The  previous  description  neglected  the  azimuthal  dependence  of  rays  arriving  and  leaving  the  foot¬ 
print.  The  full  three-dimensional  geometry  will  now  be  introduced. 

Figure  2  shows  the  geometry  of  reflection.  A  coordinate  system  was  chosen  whose  origin  is  the 
point  of  reflection  with  the  X-axis  pointing  upwind,  the  Z-axis  pointing  toward  the  zenith,  and  the 
F-axis  pointing  crosswind  such  that  a  right-handed  system  is  formed.  The  X-Y  plane  is  horizontal  at 
the  point  of  reflection.  The  tilted  facet  passes  through  the  origin. 


Z 


Figure  2.  Coordinate  system  for  facet  reflection. 

Define  a  unit  vector  U  =  (0,  (j>)  with  polar  coordinates  0,  the  zenith  angle,  and  <p,  the  azimuth.  If 
we  denote  the  Cartesian  coordinates  of  U  by  (a,  b,  c),  then  we  have 

a  =  sin  0  cos  0 

b  =  sin0sin0  (2) 

C  =  COS0 

for  the  Cartesian  coordinates  of  U  in  terms  of  its  spherical  coordinates  and 

0  =  COS“'(c) 

(j)  =  iasr\b/a) 


5 


for  the  spherical  coordinates  of  U  in  terms  of  its  Cartesian  coordinates.  Two  unit  vectors  are  shown 
in  figure  2:  U^,  pointing  from  the  origin  to  the  source,  and  U^,  pointing  from  the  origin  to  the 
receiver.  A  third  unit  vector,  U„,  is  normal  to  the  facet  at  the  point  of  reflection  but  was  removed 
from  the  figure  for  clarity^. 

The  facet  slope  in  the  upwind  direction,  tx,  is  given  by  the  slope  of  the  line  formed  at  the  intersec¬ 
tion  of  the  facet  with  the  X-Z  plane.  The  facet  slope  in  the  crosswind  direction,  ^y,  is  given  by  the 
slope  of  the  line  formed  at  the  intersection  of  the  facet  with  the  Y-Z  plane.  In  terms  of  the  Cartesian 
coordinates  of  the  facet  normal  these  slopes  are 

6.  SPECULAR 

If  a  specular  reflection  occurs,  the  three  vectors  for  source,  receiver,  and  facet  normal  are  con¬ 
nected  by  the  law  of  reflection: 

U,  4-  U,  =  2  COSO)  U„  (5) 

where  o)  is  the  angle  of  incidence  and  the  angle  of  reflection. 


bn/ 


REFLECTION 


(4) 


7.  THE  OCCURRENCE  PROBABILITY 


Following  Cox  and  Munk,  let  P  stand  for  the  probability 

P  =  piJ;xXy,W)  dtx€y 

that  a  wave  facet  will  occur  with  a  slope  within  ±  dl^x^  of  £c  ^nd  ±  d^y/2  of  ^y  when  the  wind  speed 
is  W.  The  wave  slope  occurrence  probability  density,  p,  is  proportional  to  the  horizontal  projection  of 
the  facet.  Cox  and  Munk  obtained  an  expression  for  p  whose  lowest  order  term  is 


p{^zXy,W) 


1 

InojJc 


exp 


1 

2 


al  =  0.000  -1-  3.16  •  10-^ W 
al  =  0.003  +  1.92  ■  IQ-^W 


(7) 


Here  <7^  and  <7^  variances  in  ^nd  ^y  respectively  and  IV  is  the  wind  speed  in  m  s"'.  Figure  3 

shows  the  dependence  of  p  throughout  slope  space  for  a  wind  speed  of  10  m  s'^.  The  coordinate  sys¬ 
tem  of  figure  2  has  been  inserted  at  the  top  of  the  figure  to  illustrate  the  relation  between  coordinates 
and  slopes.  Note  that  the  first  X-Y  quadrant  corresponds  to  negative  slopes. 


6.  The  zenith  angle  of  U„  is  the  same  as  the  tilt  of  the  facet.  The  tilt  is  the  angle  of  the  steepest  ascent  within  the  facet. 
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Figure  3.  Cox-Munk  occurrence  probability  density  for  a 
windspeed  of  10  m  s'^ 


8.  THE  INTERACTION  PROBABILITY 

Following  a  suggestion  of  Plass,  et  al.  (1976),  let  Q  stand  for  the  (different)  probability 

Q  =  qg;xXy,0,(t>,W)tx<^y  (8) 

that  a  facet  whose  slope  is  within  ±d^x^  of  L  and  ±dt,y/l  of  ^y  will  interact  with  a  ray  arriving 
from  the  arbitrary  direction  U  =  (6, 0)  when  the  wind  speed  is  W.  The  wave  slope  interaction  proba¬ 
bility  density,  q,  is  proportional  to  the  facet  area  projected  normal  to  the  ray.  It  has  previously  been 
shown  (Zeisse,  1994,  1995)^  that 


q{^xXy,0,4>,W) 


COSO) 

cos9„^ 


W  <  Ttfl 


COSO) 

cos6„ 


P  d^x  dty 


U  =  const. 


(9) 


Figure  4  is  a  graph  of  equation  (9),  also  for  a  wind  speed  of  10  m  s'*,  showing  how  facets  with  a  spe¬ 
cified  slope  interact  with  a  ray  pointing  in  the  direction  (80°,  270°). 


Jt 

1.  Equation  (9)  is  only  defined  for  <y  £  -j. 
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Figure  4.  Cox-Munk-Plass  interaction  probability  density  for  a 
windspeed  of  1 0  m  s'^ 

9.  EQUATIONS  FOR  SEA  RADIANCE 

The  capillary  wave  contributions  to  sea  radiance  are 


N^^ier,(l)r,W,v)  =  II  Ns{e„(j)s,v]p{(O,v)q{tx,^y,0r,4>r,W)  d^y 

0^,0)  <  Till 
Ur  ~  const. 

^ sun^ Ot  0 0'>  0  r?  ^ 

p  (a),v)  secw  sec^On  q(,txXy’^n<t>nW)  sinds  dds  d4>s 

disk 

Ur  —  const. 

Nsea(^n<PnW^Tsea,v)  =  N ^ea^v]  ||  [l  -  p[(^^v)\q{^xXy^^n^ry^)  d^x  d^y  (12) 

(0  ^  njl 
Ur  =  const. 

In  each  of  the  integrals  (10)  through  (12),  ^  plays  the  role  of  a  weighting  function  attached  to  the 
facet.  The  weight  is  applied  to  the  ray  leaving  the  footprint  and  propagating  toward  the  receiver,  and 


(10) 

(11) 
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that  ray  and  those  receiver  coordinates  are  held  constant  in  all  of  the  integrals.  A  physical  description 
and  some  mathematical  details  of  each  equation  will  now  be  presented. 


In  the  integrand  of  (10),  the  product  of  and  P  represents  the  radiance  leaving  a  single  facet 
when  Ns{6s,  <ps,  v)  is  the  spectral  sky  radiance  incident  on  that  facet  at  zenith  angle  6s  and 
azimuth  <ps-  This  product  is  weighted  by  q  and  integrated  over  all  slopes  in  the  ocean.  During  integra¬ 
tion,  a  specular  reflection  occurs  at  one  facet  after  another  inside  the  footprint  with  the  outgoing 
(receiver)  ray  held  fixed.  The  source  ray  is  swept  across  the  sky  and  sun.  Equation  (10)  will  require 
explicit  expressions  for  each  of  its  variables  in  terms  of  slopes  and  receiver  coordinates.  From  (2) 
and  (4)  it  can  be  shown  that  the  facet  tilt  is  given  in  terms  of  the  facet  slopes  by 


CQSdn  =  Cn 


1 

yi  +  e  +  ?? 


(13) 


while  the  fact  that  (o  is  the  the  angle  between  the  facet  normal  and  the  receiver  ray  implies  that 


COSO)  =  U„  •  U;. 

QfiClf  -f  b fib r  F 

~  ~  ^y6r  "I"  e^-jcn 

{-  sin^r  COS<pr  -  SinOf  sincpr  +  COS0;-} 


(14) 


Equations  (13)  and  (14)  hold  at  all  times,  regardless  of  whether  a  specular  reflection  is  taking 
place.  When  a  specular  reflection  does  occur,  the  Z  component  of  the  law  of  reflection 


U.5  =  2  costa  U„  -  U;- 


(15) 


gives 


cosds  =  2  costa  Cn  -  Cf 


2{  }  -  Cf/cl 

‘  l/c|  (>«) 

-2  sin^r  (Cx  cos<pr  -I-  sin0;-)  +  cosdf 

1  + 

where  {  }  represents  the  expression  within  curly  braces  in  (14).  Finally,  the  X  and  Y  components  of 
(15)  give 

tan0,  =  ^ 

_  2  costa  bn  -  bf 
2  COSO)  an  -  Ur 

^  2g,|  l  +  Vc;  (17) 

'^y  {  1 

(l  +  sin0,  -  {2^£y)  cos^r  +  (2^y)  cotdf 

(1  -  ^2  +  cos0r  -  (2tx^>)  sin0r  +  (2?x)  COtOf 
for  the  source  azimuth  during  specular  reflection  by  a  facet  into  a  receiver  at  (0„  4>r)- 
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Equations  (13),  (14),  (16),  and  (17)  should  be  used  in  (10)  [and  in  equation(9)  when  using  (10)]. 
The  Cartesian  expressions  are  convenient  for  computer  calculation  while  the  spherical  expressions 
are  consistent  with  the  form  of  equations  (10)  through  (12). 

In  the  integrand  of  (11),  the  product  of  Ng  and  p  represents  the  spectral  radiance  leaving  a  single 
facet  when  NgiOg,  0^,  v)  is  the  spectral  radiance  arriving  at  that  facet  from  the  sun  whose  center  is  at 
(dg,  <pg).  The  remaining  factors  in  (11),  apart  from  q,  are  the  Jacobian  of  the  transformation  from 
ocean  slopes  to  sky  coordinates  (Zeisse,  1994).  As  before,  the  integrand  is  weighted  by  q  but  now  the 
integration  is  over  the  solar  disk  in  the  sky.  (It  is  assumed  in  (1 1)  that  Ngidg,  (pg,  v)  does  not  vary  dur¬ 
ing  integration  because  the  sun  is  a  Lambertian  source  and  the  solar  disk  is  small.)  During  integra¬ 
tion,  a  specular  reflection  from  the  sun  to  the  receiver  occurs  at  those  facets  within  the  footprint  with 
the  correct  slope.  Explicit  expressions  for  each  of  the  variables  in  terms  of  source  and  receiver  coor¬ 
dinates  will  be  required  in  (11).  The  law  of  reflection 

2  coscu  U„  =  U,  -b  U,  (18) 

gives  the  facet  position  in  terms  of  the  source  and  receiver  positions  whenever  a  specular  reflection 
occurs.  The  components  of  (18)  give 


and 


I  df 

Cs  +  Cr 

sin0^  COS05  +  sin^r  cos4>r 
cosds  +  cos  dr 


(19) 


_ bs  +  br 

~  Cg  +  Cr 

sin^i  sin0^  -b  sin0^  sin^^ 

COsOs  +  cosdr 

while  its  square  gives 

2  COSTCO  =  I  +  Us  •  Ur 

=  \  +  ds  dr  +  bs  br  +  Cs  Cr 

=  1  -b  sin0i  sin0;-  cos(0j-  (pr)  +  cos05cos0r 


(20) 


(21) 


10 


Finally,  from  (2),  (4),  and  (18)  we  have 

tan^0„  =  Cx  + 

_  (Os  +  ClrY  +  {bs  +  br) 

sin20,  +  sin^6r  +  2sin05  sinS^  cos(0j  -  (pr) 

2 

(cos0^  +  COS0;-) 

Expressions  (19)  through  (22)  should  be  used  in  (11)  [and  in  (9)  when  using(ll)].  They  apply  only 
when  there  is  a  specular  reflection. 

In  (12)  there  is  no  incident  ray  or  specular  reflection,  and  integration  is  over  all  slopes  in  the 
ocean.  The  integral  in  (12)  is  the  effective  spectral  emissivity  of  the  ocean.  Explicit  expressions  in 
terms  of  slopes  and  receiver  coordinates  will  also  be  required  for  each  of  the  variables  in  (12)  [and  in 
(9)  when  using  (12)].  Equation  (13)  is  the  expression  for  and  equation  (14)  is  the  expression  for  O). 


10.  SEARAD 


SeaRad  consists  of  new  routines  added  to  MODTRAN2  to  compute  the  spectral  values  of  Nsky, 
Nsun,  and  N^ea  according  to  equations  (10),  (1 1 ),  and  (12),  respectively.  Through  modifications  to 
subroutine  “TRANS”,  these  values  are  assembled  according  to  (1)  and  integrated  over  the  wave 
number  after  obtaining  proper  path  radiance  and  transmittance  spectral  values.  SeaRad  also 
introduces  minor  changes  in  subroutine  “DPFNMN”  and  major  changes  in  subroutine  “DRIVER  . 
These  changes  will  now  be  considered  in  more  detail. 

The  modifications  to  “DRIVER”  are  briefly  shown  in  figure  5.  After  the  normal  call  to  “GEO”,  a 
test  is  conducted  to  see  whether  the  ray  chosen  by  the  user  has  hit  the  surface  of  the  sea.  If  so,  geom¬ 
etry  cards  required  for  the  sea  calculation  are  issued  by  subroutine  “Card”  to  file  “TapeS.Sea”,  and 
input  is  temporarily  redirected  to  “TapeS.Sea”.  An  example  of  “TapeS.Sea”  is  given  on  page  A-4,  for 
the  run  initiated  by  the  input  file  on  page  A-2.  After  the  final  card  has  been  read  fi-om  “TapeS.Sea”, 
sea  radiance  is  calculated  in  “TRANS”.  Then  “TAPES”  is  restored  as  the  active  input  file  and  normal 
program  execution  is  resumed.  Please  see  Appendix  B  for  a  detailed  flowchart  of  the  modifications 
to  “DRIVER”  as  well  as  the  complete  source  code  for  the  modified  version  of  “DRIVER”. 

Conditions  in  “DPFNMN”  determine  whether  or  not  the  sea  has  been  hit.  “DPFNMN”  is  a  subrou¬ 
tine  reached  by  a  sequence  of  calls  beginning  in  the  driver  with  a  call  to  subroutine  “GEO”.  Modifi¬ 
cations  to  “DPFNMN”  are  summarized  in  figure  6.  A  logical  variable  “Sea”,  initially  set  false,  is  set 
true  in  “DPFNMN”  if  the  following  four  conditions  are  met: 

1 .  The  program  has  reached  the  section  of  code  following  the  comment  line  “Tangent  path  inter¬ 
sects  earth.” 

2.  The  user  has  chosen  a  radiance  mode. 

3.  The  variable  “HMIN”  is  equal  to  zero. 

4.  The  variable  “SeaSwitch”  is  true. 
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Figure  5.  Flowchart  for  modifications  to  MODTRAN2 
subroutine  “DRIVER.” 
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Figure  6.  Flowchart  for  modifications 
to  MODTRAN2  subroutine  “DPFNMN.” 

The  variable  “Sea”  is  stored  in  a  common  block  made  available  to  the  driver,  which  inspects  “Sea” 
before  and  after  each  of  its  calls  to  “GEO”.  A  change  from  false  to  true  indicates  that  the  ocean  has 
been  hit  during  that  call.  A  hit  induces  a  geometry  calculation  by  a  call  to  subroutine  “Foot  (if 
EEMSCT  =  1)  or  subroutine  “SunFoot”  (if  ffiMSCT  =  2).  This  is  followed  in  each  case  by  a  call  to 
subroutine  “Card”. 

The  purpose  of  “Card”  is  to  supply  sources  for  the  Cox-Munk  routines  “Sky”  and  “Sun.”  As 
shown  in  figure  7,  geometry  cards  are  issued  here  to  file  “TapeS.Sea”  to  obtain  spectral  radiance 
along  paths  to  the  sky  at  three  separate  zenith  angles.  These  three  cards,  one  for  each  zenith  angle, 
are  called  “Sky  Cards”  in  the  flowchart.  Later  these  data  will  be  used  by  subroutine  “Fit”  to  establish 
a  two-parameter  least  squares  fit  at  each  wave  number  providing  “Sky”  with  the  sky  dome  radiance. 
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Figure?.  Flowchart  for 
SeaRad  subroutine  “Card.” 
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If  the  sun  is  involved,  “Card”  will  issue  a  fourth  and  final  Card  3,  called  a  Sun  Card  in  the  flow¬ 
chart,  which  provides  solar  irradiance  for  later  use  as  a  source  by  subroutine  “Sun”. 

The  modifications  described  to  this  point  have,  in  effect,  inserted  three  sky  cards  (followed  by  a 
sun  card  if  necessary)  into  the  user’s  input  file  without  the  user’s  knowledge.  The  insertion  is  made 
only  if  the  user  has  chosen  a  Card  3  whose  path  terminates  on  the  surface  of  the  earth.  Such  a  Card  3 
is  called  a  “Path  Card”  in  figure  C-1.  Within  the  wave  number  integration  loop  in  “TRANS”,  spectral 
values  of  transmission,  incident  sky  radiance,  and  incident  solar  irradiance  are  stored  in  arrays 
Tau(V),  Nsky(V),  and  Ho(V),  respectively.  Outside  the  wave  number  integration  loop  these  values 
are  recalled  for  the  sea  radiance  calculation  by  subroutine  “Sky”  (or  subroutine  “Sun”  if  lEMSCT 
=  2). 

The  modified  version  of  “DRIVER”  is  contained  in  Appendix  B  along  with  a  detailed  flowchart  of 
its  modifications.  Appendix  C  contains  the  source  code  and  a  flowchart  for  the  modified  version  of 
“TRANS”,  and  Appendix  D  contains  new  code  for  the  sea  radiance  calculation. 

11.  CONCLUSION 

SeaRad,  a  modification  of  MODTRAN2,  computes  sea  radiance  between  52.63  cm"’  and 
25000  cm-f  Preliminary  comparisons  with  data  show  that  SeaRad  has  an  approximate  accuracy  of 
several  °C  in  the  infrared. 

SeaRad  is  currently  designed  for  a  single  pixel  and  takes  approximately  10  s  to  execute.  Each  time 
a  new  geometry  is  chosen  by  the  user,  SeaRad  recalculates  the  source  radiance  and  the  path  radiance 
and  transmission.  However,  only  the  path  properties  change  significantly  from  one  pixel  to  the  next 
in  a  typical  ocean  image.  If  SeaRad  were  redesigned  to  apply  to  sea  images,  the  speed  per  pixel  could 
be  reduced,  up  to  a  factor  of  almost  four,  by  calculating  values  of  source  radiance  just  once  for  the 
entire  image. 


15 


12.  REFERENCES 


Berk,  A.,  L.  S.  Bernstein,  and  D.  C.  Robertson.  1989.  “MODTRAN;  A  Moderate  Resolution  Model 
for  LOWTRAN  7,”  Report  GL-TR-89-0122.  Air  Force  Geophysics  Laboratory,  Bedford,  MA. 

Cox,  C.  and  W.  Munk.  1954.  “Measurement  of  the  Roughness  of  the  Sea  Surface  from  Photographs 
of  the  Sun’s  Glitter,”  Journal  of  the  Optical  Society  of  America  44:  838-850. 

Cox,  C.  and  W.  Munk.  1956.  “Slopes  of  the  Sea  Surface  Deduced  from  Photographs  of  Sun  Glitter,” 
Scripps  Institution  of  Oceanography  Bulletin  6:  401-487. 

Fenn,  R.  W.,  S.  A.  Clough,  W.  O.  Gallery,  R.  E.  Good,  F.  X.  Kneizys,  J.  D.  Mill,  L.  S.  Rothman, 

E.  P.  Shettle,  and  F.  E.  Volz.  1985.  “Optical  and  Infrared  Properties  of  the  Atmosphere,”  In 
Handbook  of  Geophysics  and  the  Space  Environment.  A.  S.  Jursa,  Ed.  Air  Force  Geophysics 
Laboratory,  Bedford,  MA. 

Hale,  G.  M.  and  M.  R.  Querry.  1973.  “Optical  Constants  of  Water  in  the  200  nm  to  200  [tm  Wave¬ 
length  Region,”  Applied  Optics  3:  555-563. 

Isaacs,  R.  G.,  W.  C.  Wang,  R.  D.  Worsham,  and  S.  Goldenberg.  1986.  “Multiple  Scattering  Treatment 
for  Use  in  the  LOWTRAN  and  FASCODE  Models,”  Report  AFGL-TR-86-0073:  Air  Force 
Geophysics  Laboratory,  Bedford,  MA. 

Kneizys,  F.  X.,  E.  P.  Shettle,  L.  W.  Abreu,  J.  H.  Chetwynd,  G.  P.  Anderson,  W.  O.  Gallery, 

J.  E.  A.  Selby,  and  S.  A.  Clough.  1988.  “Users  Guide  to  LOWTRAN  7,”  Report  AFGL- 
TR-88-0177.  Air  Force  Geophysics  Laboratory,  Bedford,  MA. 

Kneizys,  F.  X.,  E.  P.  Shettle,  W.  O.  Gallery,  J.  H.  Chetwynd,  L.  W.  Abreu,  J.  E.  A.  Selby, 

S.  A.  Clough,  and  R.  W.  Fenn.  1983.  “Atmospheric  Transmittance/Radiance:  Computer  Code 
LOWTRAN  6,”  Report  AFGL-TR-83-0187.  Air  Force  Geophysics  Laboratory,  Bedford,  MA. 

Kneizys,  F.  X.,  E.  P.  Shettle,  W.  O.  Gallery,  J.  H.  Chetwynd,  L.  W.  Abreu,  J.  E.  A.  Selby,  R.  W.  Fenn, 
and  R.  A.  McClatchey.  1980.  “Atmospheric  Transmittance/Radiance:  Computer  Code  LOW¬ 
TRAN  5,”  Report  AFGL-TR-80-0067.  Air  Force  Geophysics  Laboratory,  Bedford,  MA. 

Kneizys,  F.  X.,  J.  E.  A.  Selby,  J.  H.  Chetwynd,  and  R.  A.  McClatchey.  1978.  “Atmospheric  Transmit¬ 
tance/Radiance:  Computer  Code  LOWTRAN  4,”  Report  AFGL— TR— 78— 0053.  Air  Force  Geo¬ 
physics  Laboratory,  Bedford,  MA. 

Plass,  G.  N.,  G.  W.  Kattawar,  and  J.  A.  Guinn.  1975.  “Radiative  Transfer  in  the  Earth’s  Atmosphere 
and  Ocean:  Influence  of  Ocean  Waves  f  Applied  Optics  14:  1924—1936. 

Querry,  M.  R.,  W.  E.  Holland,  R.  C.  Waring,  L.  M.  Earls,  and  M.  D.  Querry.  1977.  “Relative  Reflec¬ 
tance  and  Complex  Refractive  Index  in  the  Infrared  for  Saline  Environmental  Waters,  ”  Journal 
of  Geophysical  Research  82:  1425-1433. 

Selby,  J.  E.  A.,  E.  P.  Shettle,  and  R.  A.  McClatchey.  1976.  “Atmospheric  Transmittance  from  0.25  to 
28.5  [xm:  Supplement  LOWTRAN  3B,”  Report  AFGL-TR-76-0258.  Air  Force  Geophysics  Lab¬ 
oratory,  Bedford,  MA. 

Selby,  J.  E.  A.  and  R.  A.  McClatchey.  1975.  “Atmospheric  Transmittance  from  0.25  to  28.5  ftm: 
Computer  Code  LOWTRAN3,”  Report  AFGL-TR-72-0255.  Air  Force  Geophysics  Laboratory, 
Bedford,  MA. 


16 


Selby,  J.  E.  A.  and  R.  A.  McClatchey.  1972.  “Atmospheric  Transmittance  form  0.25  to  28.5  [xm: 
Computer  Code  LOWTRAN  2,”  Report  AFGL-TR-72-0745.  Air  Force  Geophysics  Laboratory, 
Bedford  MA. 

Stratton,  J.  A.  1941.  Electromagnetic  Theory,  pp.  505  ff.,  McGraw-Hill,  New  York,  NY. 

Zeisse,  C.  R.  1994.  “Radiance  of  the  Ocean  Horizon,”  NRaD  TR  1660  (April).  Naval  Command, 
Control  and  Ocean  Surveillance  Center,  RDT&E  Division,  San  Diego,  CA. 

Zeisse,  C.  R.  1995.  “Radiance  of  Ocean  Horizon,”  Journal  of  the  Optical  Society,  vol.  A,  pp. 
2022-2030. 


17 


APPENDIX  A 

Sea/?ac/ INPUT  AND  OUTPUT  FILES 


A-1 


w  o\ 


“TAPESRAD.STD”  INPUT  FILE 


3 
0 

00.023 

940 

0 


C;\MOD2\TAPE5RAD.FIL  7/20/95 


1  1 

0  5 

.  000 
950 


0  0 

0  0 

100.000 

10 


0  0 
10.000 
.000 
5 


0  0 
10.000 
.  000 


0  0 
10.000 
0.00 


0  288.15  0.00 

.000  ,000 
0  90.000  T 


A-2 


“OUT”  FILE 


C: \MOD2\OUT,FIL  7/20/95 

*****  SEAIUVD,  A  MODIFICATION  OF  LOWTRAN7  ***** 

DATE:  07/20/95  TIME:  15:11:38 

THERMAL  RADIANCE  MODE 

MULTIPLE  SCATTERING  USED 

MARINE  AEROSOL  MODEL  USED 

WIND  SPEED 
WIND  SPEED 
RELATIVE  HUMIDITY 
AIRMASS  CHARACTER 
VISIBILITY 

SLANT  PATH  TO  SPACE 

HI  =  0.023  KM 

HMIN  =  0.000  KM 

ANGLE  =  100.000  DEG 

FREQUENCY  RANGE 

IVl  =  940  CM->1 

IV2  =  950  CM-l 

IDV  =  10  CM-1 

IFWHM  =  5  CM-1 

IFILTER  =  0 

SUMMARY  OF  THE  GEOMETRY  CALCULATION 


HI 

= 

0.023 

KM 

H2 

= 

0.000 

KM 

ANGLE 

= 

100.000 

DEG 

RANGE 

0.133 

KM 

BETA 

0.001 

DEG 

PHI 

= 

80.001 

DEG 

HMIN 

0.000 

KM 

BENDING 

= 

0.000 

DEG 

LEN 

= 

0 

SEA  AT  288.15  K  REPLACES  BLACK  BODY  BOUNDARY 

UPWIND  =  90.000  DEG  EAST  OF  LINE  OF  SIGHT 


(  10.64  MI CROMETERS ) 

(  1 0 . 5  3  MI CROMETERS ) 


10.00  M/SEC 

10.00  M/SEC,  24  HR  AVERAGE 
50.00  PERCENT 
5 

10.00  KM 


RECEIVED  RADIANCE  VALUES 


PATH  TO  FOOTPRINT  0 
SEA  EMISSION  =  0 
SKY  REFLECTION  =  0 
SUN  GLINT  =  0 

TOTAL  RADIANCE  =  0 
BLACK  BODY  TEMP. 


01814  W  M-2  SR-1  (AV.  TRANS.  0.9776) 
70712  W  M-2  SR-1 
06125  W  M-2  SR-1 
00000  W  M-2  SR-1 

78652  W  M-2  SR-1 
6.7  C 


A"3 


“TAPE5.SEA”  FILE 


C:\MOD2\TAPE5SEA.FIL  7/20/95 


0.000 

0.000 

57.296 

0.000 

0.000 

0.000 

0 

90.000 

T 

0.000 

0.000 

73.148 

0.000 

0.000 

0.000 

0 

90.000 

T 

0.000 

0.000 

89.000 

0.000 

0.000 

0.000 

0 

90.000 

T 

A-4 


APPENDIX  B 

MODIFIED  SUBROUTINE  “DRIVER” 


B-1 


BEGIN 

MODTRAN 


INITIALIZE 


SEA=.F. 

SEAOLD=.F. 

HIT=.F. 

DONE=.F. 

MSEA=-1 


READ 
CARD  1 


SAVE  CARD  1 
PARAMETERS 


*old=* 

(*=ITYPE. 

lEMSCT 

IMULT, 

TBOUND, 

SALB.) 


READ 
CARD  2 


CONTINUE 


^  FIRST  SKY  ^ 
^CARD  NEXT?^ 


ITYPE=3 

SALB=0 


SUN  CARD 
.  NEXT?  . 


IEMSCT=3 

IMULT=0 

LSAME=,F. 

LFIRST=.T. 


Figure  B-1.  Detailed  flowchart  for  modified  subroutine  “DRIVER”. 
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Figure  B-1 .  Detailed  flowchart  for  modified  subroutine  “DRIVER”.  (Continued) 
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400 


Figure  B-1 .  Detailed  flowchart  for  modified  subroutine  “DRIVER”.  (Continued) 


Figure  B-1 .  Detailed  flowchart  for  modified  subroutine  “DRIVER”.  (Continued) 


B-5 


c 


c 


SUBROUTINE  DRIVER 

COMMON  RELHUM(34)  ,HSTOR(34)  ,ICH(4)  ,VH(17)  ,TX(63)  ,W(63) 

COMMON  WRATH ( 68 ,  63 ) , TBBY (68 ) 

COMMON  IMSMX, WRATH (102, 63)  , TBBY (102)  ,RATM(102)  ,NSPEC,KPOINT (12) 
COMMON  ABSC(5,47)  ,EXTC(5,47) ,ASYM(5,47)  ,VX2(47) ,AWCCON(5) 

COMMON  /IFIL/  IRD,IPR,IPU,NPR,IPR1,IP6,IP7,IP8,IR4,IRDS,IP6S,ITR, 
+  Isky,Isun,Ipath 

COMMON  /CARDl/  MODEL,  ITYPE,  IEMSCT,M1,M2  ,M3  ,  IM,NOPRT,TBOUND,  SALB 
1  , MODTRN 

LOGICAL  MODTRN 
logical  ground 
logical  Isaiae 

LOGICAL  SeaSwitch , Sea , SeaOld , Hit , Done 
COMMON  /CARDIA/  M4  ,M5  ,M6,MDEF,  IRDl ,  IRD2 

COMMON  /CARD2/  IHAZE,  ISEASN,  IVULCN ,  ICSTL,  ICLD ,  IVSA,  VIS ,  WSS  ,  WHH, 

1  RAINRT 

COMMON  /CARD2A/  CTHIK , CALT , CEXT 
COMMON  /CARD2D/  IREG (4 ) , ALTB (4 ) , IREGC (4 ) 

COMMON  /CARD3/  HI  ,H2  ,  ANGLE, RANGE,  BETA, RE,  LEN,Psi ,  SeaSwitch 
COMMON  /CardSAl/  IPARM,IPH,IDAY,ISOURC 

COMMON  /Card3A2/  PARMl ,  PARM2  ,  PARM3  , PARM4  , GMT, PSIPO,  ANGLEM,  G 
COMMON  /CARD4/  IVl , IV2 , IDV, IFWHM, IFILTER 
COMMON  /CNSTNS/  PIX,  CA,DEG,GCAIR,BIGNUM,BIGEXP 
COMMON  /CNTRL/  KMAX,M,  IKMAX,  NL,  ML,  IKLO,  ISSGEO ,  IMULT 
COMMON  /MODEL/  ZM(34) , PM(34) ,TM(34) ,RFNDX (34) , DENSTY (63 , 34) , 

1  CLDAMT(34)  ,RRAMT(34)  ,EQLWC(34)  ,HAZEC(34) 

COMMON  /SOLS/  AHl (68)  , ARH(68)  , 

1  WPATHS(102, 63)  ,PA(68)  ,PRX(68)  ,ATHETA(35)  ,ADBETA(35)  ,LJ(69)  , 

2  JTURN,ANGSUN,CSZEN(68)  ,TBBYS  (102 , 12)  ,  PATMS  (102 , 12 ) 

COMMON  /MART/  RHH 

COMMON  /USRDTA/  NANGLS,ANGF(50) , F (4 , 50) 

COMMON  /MDLZ/  HMDLZ(8) 

COMMON  /ZVSALY/  ZVSA(IO)  ,RHVSA(10)  ,  AHVSA(IO)  ,  IHVSA(IO) 

CHARACTER*4  HHAZE ,  HSEASN , HVULCN ,  BLANK , HMET ,  HMODEL ,  HTRRAD 
COMMON  /TITL/  HHAZE(5, 16) ,HSEASN(5, 2)  ,HVULCN(5, 8) ,BLANK, 

1  HMET(5,2) ,HMODEL(5,8) ,HTRRAD(6,4) 

COMMON  /VSBD/  VSB(IO) 

COMMON  /MNLT/TBBSS(68) ,TBBMS(34) ,WPMS(34,63) , IMSMX, WPMSS (34 , 63 ) 
COMMON  /PATH/PL (68) ,QTHETA(68) , ITEST, HI , HF, AHT ( 68 ) ,tph(68) 

COMMON  /AERTM/TAE7  ,  TAE12  ,  TAE13  ,  TAE14  ,  TAE16 

common  /graund/gndalt 

common  /small3/ small 

common  /solar/ Is ame 

PARAMETER  (Kr  =  216,  Kv  ==  400) 

COMMON  /Filters/  FLIST(5,6), 

+  F1LTER1(45),  BBl(Kr),  FILTER2(54),  BB2 (Kr) , 

+  FILTER3(39),  BB3 (Kr) ,  FILTER4(47),  BB4 (Kr) , 

+  FILTERS ( 101 ), BBS (Kr) ,  FILTER6(75),  BB6(Kr) 

COMMON  /Constants/  pi, r2d,d2r, epsilon, delta, onem, onep, inf inity 
COMMON  /Geometry/  Tsun, Psun, Tr , Pr 
COMMON  /Sea/  Sea, Hit,Msea,TBOUNDold, lEMSCTold 
COMMON  /Sealndex/  AlphaOl (100) ,  Alpha02(20), 

+  BetaOl  UOO)  ,  Beta02  (20) 

logical  Ifirst 
data  If irst/ . true. / 

REAL  infinity 
CHARACTER*8  Date$,  Time$ 
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CHARACTER* 14  Prog$ 

INTEGER*4  Istart,  lend 

C  First,  get  starting  time  to  measure  total  execution  time: 

CALL  TIMER (Istart) 
c 

c  Ifirst  is  true  when  first  solar  parameters  are  read  in  a  series 

c  of  runs  involving  solar  parameters, 

c 

C*****HDATE  AND  HTIME  CARRY  THE  DATA  AND  TIME  AND  MUST  BE  DOUBLE 
C*****PRECISION  ON  A  32  BIT  WORD  COMPUTER 
C@  DOUBLE  PRECISION  HDATE, HTIME 

DIMENSION  PLST(68) ,CSENSV(68) ,QTHETS(68) 

DATA  IRPT  /  0  / 

C*****IRD,  IPR,  AND  IPU  ARE  UNIT  NUMBERS  FOR  INPUT,  OUTPUT,  AND 
C*****IPR1  =  OUTPUT  OF  MOLECULAR  TRANSMITTANCE 
DATA  MAXGEO  /  68/ 

small  =  2.0 

IP4  =  14 
IRD  =  5 
IPR  =  6 
IP6  =  16 
IP6S=  26 
IPU  =  7 
IP7  =  17 
IPR1=  8 
IPS  =  18 
IRDS  “29 
ITR  =  30 
ISCRCH  =  10 
ITM  =  31 
Isky  =  32 
Isun  =33 
Ipath  =34 

OPEN  (IRD,  FILE= *  TAPES  * ,  STATUS= * OLD  * ) 

OPEN  (IPR,  FILE=*TAPE6* ,  STATUS= * UNKNOWN ' ) 

OPEN  ( IP6 ,  FILE=  *  OUT  * ,  STATUS=  *  UNKNOWN  » ) 

OPEN  ( IPU ,  FILE=  *  TAPE7  » ,  STATUS=  *  UNKNOWN  * ) 

OPEN  ( IP7 ,  FILE= • TAPE7 . PLT  * , STATUS= » UNKNOWN  * ) 

OPEN  ( IPRl , FILE=  » TAPES  » ,  STATUS=  *  UNKNOWN  * ) 

OPEN  (IPS,  FILE= 'TAPES. PLT * ,STATUS=* UNKNOWN*) 

OPEN  ( IP4 ,  FILE=  *  OUT , PLT  * ,  STATUS=  *  UNKNOWN  * ) 

OPEN  ( ITR,  FILE=  *  TRANS . PLT  * , STATUS= ' UNKNOWN  * ) 

OPEN  ( ISCRCH , STATUS=  *  SCRATCH  * , FORM=  *  UNFORMATTED » ) 

OPEN  (Isky,  FILE=*Sky.plt* ,  STATUS= » UNKNOWN  * ) 

OPEN  ( I sun ,  FILE=  *  Sun . pit  * ,  STATUS=  *  UNKNOWN  * ) 

OPEN  ( Ipath , FILE=  *  Path . pit  * , STATUS=  *  UNKNOWN  * ) 

OPEN  ( ITM,  FILE=  *  TIME  » ,  STATUS=  *  UNKNOWN  * ) 

C 

C  ALTITUDE  PARAMETERS 

C 

C  ZMDL  COMMON/MODEL/  THE  ALTITUDES  USED  IN  LOWTRAN 

C  ZCVSA,ZTVSA,ZIVSA  CARD  3.3  LOWTRAN  FOR  VSA  INPUT 

C  ZVSA  NINE  ALTITUDES  GEN  BY  VSA  ROUTINE 

C 

Pix=2.0*ASIN(1.0) 

CA=Pix/180. 
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driv  680 
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driv  700 
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DEG=  1.0/CA 


pi 

r2d 

d2r 

epsilon 

delta 


=  Pix 
=  180-/pi 
=  pi/180. 

=  d2r*0.2659 
=  1.4E-6 
=  1.  ^  delta 
=  1.  +  delta 
=  999999 


onep  =1.  +  aeira 

infinity  =  999999 

RANGE=0.0  dr IV  870 

C*****gCAIR  is  the  gas  constant  for  air  in  units  of  MB/(GM  CM-3  K)  driv  880 

GCAIR  =  2.87053E+3  890 

C*****BIGNUM  AND  BIGEXP  ARE  THE  LARGEST  NUMBER  AND  THE  LARGEST  ARGUMENT  driv  900 
C*****EXP  ALLOWED  AND  ARE  MACHINE  DEPENDENT.  THE  NUMBERS  USED  HERE  ARE  Fdriv  910 


C*****A  TYPICAL  32  BIT-WORD  COMPUTER.  920 

BIGNUM  =  1.0E35  ^^IV  930 

BIGEXP  =87.0  ^40 

C  THE  VALUES  FOR  BIGNUM  AND  BIGEXP  FOLLOW  THE  driv  950 

C  DESCRIPTION  UNDER  EXP  FUNCTION  IN  "IBM  SYSTEM  360/  driv  960 

C  AND  SYSTEM  370  FORTRAN  IV  LANGUAGE"  driv  970 

C  BIGNUM  =  4.3E68  driv  980 

C  BIGEXP  =  174.6  driv  990 

KMAX=63  drivlOOO 

C*****NL  IS  THE  NUMBER  OF  BOUNDARIES  IN  THE  STANDARD  MODELS  1  TO  6  drivlOlO 

C*****BOUNDARY  34  (AT  99999  KM)  IS  NO  LONGER  USED  drivl020 

NL  =  33  dnvl030 

********************************************************************* 

Sea  =  . FALSE - 

SeaOld  =  .FALSE. 

Hit  =  .FALSE. 

Msea  =  -1 

Done  =  .FALSE. 

********************************************************************* 

C*****CALL  TIME  AND  DATE:  drivl040 

C*****THE  USER  MAY  WISH  TO  INCLUDE  SUBROUTINES  FDATE  AND  FCLOCK  WHICH  drivlOSO 
C*****RETURN  THE  DATE  AND  TIME  IN  MM/DD/YY  AND  HH.MM.SS  FORMATS  drivl060 

C*****RESPECTIVELY.  THE  REQUIRED  ROUTINES  FOR  A  CDC  6600  ARE  INCLUDED  ATdrivl070 
C*****THE  main  PROGRAM  IN  COMMENT  CARDS.  drivlOSO 

C@  CALL  FDATE (HDATE)  drivl090 

ra  r-aT.T.  ■prr.nrx  rwTTM'R^  drivllOO 


drivl040 

drivlOSO 

drivl060 


C@  CALL  FDATE (HDATE) 

C@  CALL  FCLOCK (HTIME) 

CALL  DATE  (Date$) 

CALL  TIME  (Tiine$) 

C 

C*****START  CALCULATION 

C 

C 

100  DO  10  II  =  1,4 

10  IREG(II)  =  0 

WRITE(IPR, 1000) 

1000  FORMAT( *  1  * ,20X,  *  *****  MODTRAN  ******) 

C§  WRITE(IPR,1010)  HDATE, HTIME 

1010  FORMAT ( ’ 1  * , 20X, *  *****  MODTRAN  ***** ' , lOX, 2 (IX, A8, IX) ) 
DO  80  1=1,4 

DO  80  J=l,40 

ABSC(I, J)=0. 

EXTC(I, J)=0. 

80  ASYM(I,J)=0. 
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JPRT  =  0  drivl260 

IKL0=1  drivl270 

Q  drivl280 

C*****CARD  1  drivl290 

C  drivl300 

READ(IRD,  * (L1,I4,12I5,F8,3,F7.2)  *) MODTRN, MODEL, IT YPE, 

+  I EMSCT ,  IMULT ,  Ml ,  M2 ,  M3  ,  M4 ,  M5 ,  M6 ,  MDEF ,  IM ,  NOPRT ,  TBOUND ,  SALE  dr  ivl 340 

1110  FORMAT(13I5,F8,3,F7.2)  drivl350 


**********************************  Save  parameters  to  restore  them 
ITYPEold  =  ITYPE 
lEMSCTold  =  lEMSCT 
IMULTold  =  IMULT 
If  (TBOUND  .EQ.  0.)  TBOUND  =0.1 
TBOUNDold  =  TBOUND 
SALBold  =  SALE 

***********************************  in  case  new  geometry  cards  are 
C  later  introduced  via  file  TAPE5.SEA 


1018 

1020 


IF  (MODTRN)  THEN 

Prog$  =  »MODTRAN2  ****** 

ELSE 

Prog$  =  'LOWTRAN7  ****** 

END  IF 

WRITE  {IP6,  1018)  Prog$ 

FORMAT(15X,  ******  SEARAD,  A  MODIFICATION  OF  *,  A14) 

WRITE  (IP6,  1020)  Date$,  Time$ 

FORMAT  (/,  »DATE:*,  IX,  A8,  T60,  'TIME:*,  IX,  A8) 

SELECT  CASE  (lEMSCT) 

CASE  (0) 

WRITE  (IP6,  *(/,  18HTRANSMITTANCE  MODE)  * ) 

CASE  (1) 

WRITE  (IP6,  *(/,  2 IHTHERMAL  RADIANCE  MODE) *) 

CASE  (2) 

WRITE  (IP6,  *(/,  32HTHERMAL  PLUS  SOLAR  RADIANCE  MODE)*) 

CASE  ( 3 ) 

WRITE  (IP6,  •(/,  21HSOLAR  IRRADIANCE  MODE)  *) 

END  SELECT 


SELECT  CASE  (IMULT) 

CASE  (0) 

PRINT  *,  "IMULT  =  ",  IMULT,  ":  BEWARE  OF  BEN-SHALOM" 
WRITE  (IP6,  *(/,  22HSINGLE  SCATTERING  USED)*) 

CASE  (1) 

WRITE  (IP6,  •(/,  24HMULTIPLE  SCATTERING  USED)*) 

END  SELECT 


C 


1111 

C 

c 

c 

c 

c 

c 


WRITE  (IPR,  *  (15H0  CARD  1  *****,  LI ,  14 , 1215 ,  F8 . 3  ,  F7 . 2 )')  MODTRN, MODELdrivl380 

1  ,  ITYPE ,  lEMSCT ,  IMULT ,  Ml ,  M2 ,  M3  , M4  ,  M5 ,  M6 , MDEF ,  IM ,  NOPRT ,  TBOUND ,  SALB  dr ivl3 9 0 


FORMAT(*0  CARD  1  ******, 1315 , F8 . 3 , F7 . 2 )  drivl400 

IF  (IMULT  -EQ.  1  .AND.  NOPRT.  EQ.  1)  NOPRT  =  0  drivl410 

drivl420 

SET  THE  NUMBER  OF  SPECIES  TREATED  WITH  THE  1  CM“1  BAND  MODEL.  drivl430 

ALSO,  FOR  EACH  SPECIES,  SET  THE  POINTER  WHICH  MAPS  THE  HITRAN  drivl440 

NUMERICAL  LABEL  TO  THE  LOWTRAN  NUMERICAL  LABEL.  drivl450 

drivl460 


NSPEC=12 
KPOINT(  1)=17 
KPOINT(  2) =3 6 


drivl470 

drivl480 

drivl490 


KPOINT(  3) =31 
KPOINT(  4) =47 
KPOINT(  5) =44 
KPOINT(  6) =46 
KPOINT(  7) =50 
KPOINT(  8) =54 
KPOINT(  9) =56 
KPOINT(10)=55 
KPOINT(ll)=52 
KPOINT(12)=ll 
C 

IRDl  =  0 
IRD2  =  0 

IF  (MODEL. EQ.O)  LEW  =  0 

IF( (MODEL. EQ.O)  .OR.  (MODEL. EQ. 7) )  GO  TO  110 
IF (Ml. EQ.O)  Ml=MODEL 
IF (M2. EQ.O)  M2=MODEL 
IF (M3. EQ.O)  M3=MODEL 
IF (M4. EQ.O)  M4=MODEL 
IF (M5. EQ.O)  M5=M0DEL 
IF (M6. EQ.O)  M6=M0DEL 
IF(MDEF.EQ.O)  MDEF=1 
110  CONTINUE 
M=MODEL 
NPR  =  NOPRT 

C*****CARD  2  AEROSOL  MODEL 

READ  ( IRD ,  12  00 )  IHAZE ,  I SEASN ,  IVULCN ,  ICSTL ,  ICLD ,  IVSA ,  VIS ,  WSS ,  WHH , 

1  RAINRT,GNDALT 

1200  FORMAT(6I5,5F10.3) 

WRITE  ( IPR ,  12 01 )  IHAZE ,  ISEASN ,  IVULCN ,  ICSTL ,  ICLD IVSA ,  VIS ,  WSS ,  WHH , 
1  RAINRT,GNDALT 

IF(GNDALT.GT.O. )  WRITE(IPR, 1199) GNDALT 
1199  FORMAT (IHO,'  GNDALT  =', FIO. 2) 

IF (GNDALT. GE. 6.0)  THEN 

WRITE ( IPR , 12  0  2 ) GNDALT 
GNDALT=0. 

ENDIF 

1201  F0RMAT('0  CARD  2  *****', 615, 5F10 . 3) 

1202  FORMAT('0  GNDALT  GT  6.0  RESET  TO  ZERO,  GNDALT  WAS', FIO. 3) 

C 

IF(VIS.LE. 0.0. AND. IHAZE. GT.O)  VIS=VSB (IHAZE) 

RHH=  0. 

IF(MODEL.EQ.O.OR.MODEL.EQ.7)  GO  TO  205 

IF( (MODEL. EQ. 3. OR. MODEL. EQ. 5) .AND. ISEASN. EQ.O)  ISEASN=2 
C 

IF(IVSA.EQ.l  .AND.  IHAZE. EQ. 3) 

1  CALL  MARINE(VIS, MODEL, WSS, WHH, ICSTL, EXTC,ABSC,1) 

ICH(1)=IHAZE 

ICH(2)=6 

ICH(3)=9+IVULCN 

205  IF(RAINRT.EQ.O)  GO  TO  210 
WRITE(IPR,1205)  RAINRT 

1205  FORMATCO  RAIN  MODEL  CALLED,  RAIN  RATE  =  ',F9.2,'  MM/HR') 

210  ICH(4)=18 

IF(ICH(1) .LE.0)ICH(1)=1 
IF(ICH(3) .LE.9)ICH(3)=10 
IF(ICLD.GE.l  .AND.  ICLD.LE.il)  THEN 
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ICH(4)=ICH(3) 

ICH(3)=ICH(2) 

ICH(2)=ICLD 
END  IF 
IFLGA=0 
IFLGT=0 
CTHIK=-99. 

CALT=*-99. 

CEXT=-99. 

ISEED=-99 

IF(ICLD  .LT.  18)  GO  TO  230 
C*****CARD  2A  CIRRUS  CLOUDS 

READ ( IRD ,1210) CTHIK , CALX , CEXT , ISEED 

1210  FORMAT(3F10.3,I10) 

WRITE ( I PR , 1 2 1 1 ) CTHIK , CALT , CEXT , I S EED 

1211  FORMAT(*0  CARD  2A  **★**', 3F10 • 3 , 110) 

230  CONTINUE 

C*****CARD  2B  VERTICAL  STRUCTURE  ALGORITHM 
ZCVSA=-99. 

ZTVSA=-*99 . 

ZINVSA=-99. 

C 

IF(  IVSA.  EQ.  0  )  GO  TO  240 
READ  (IRD, 1230)  ZCVSA, ZTVSA, ZINVSA 

1230  FORMAT  (3F10. 3) 

WRITE  ( IPR  ,1231)  ZCVSA ,  ZTVSA ,  ZINVSA 

1231  FORMAT(*0  CARD  2B  ******, 3F10 . 3 ) 

C 

CALL  VSA  ( IHAZE,  VIS ,  ZCVSA,  ZTVSA,  ZINVSA,  ZVSA , RHVSA ,  AHVS A ,  IHVSA) 
C 

C  END  OF  VSA  MODEL  SET-UP 

C 

240  IF  (MODEL. NE.O  .AND.  MODEL. NE. 7  )  ML=NL 
MDELS=MODEL 
DO  250  1=1,5 

IF (MDELS . NE . 0 ) HMODEL (1,7) =HMODEL ( I , MDELS ) 

250  IF (MDELS . EQ . 0 ) HMODEL (1,7) =HMODEL (1,8) 

C 

IF(IM  .EQ.  1)  THEN 

IF( (MODEL.EQ.7.AND.IM.EQ.1)  .OR. (MODEL. EQ.O) )  THEN 
C 

C*****CARD  2C  USER  SUPPLIED  ATMOSPHERIC  PROFILE 
C 

READ  (IRD, 1250)  ML, IRDl, IRD2 , (HMODEL( I , 7 ) , 1=1 , 5) 

1250  FORMAT (315, 18A4) 

WRITE ( IPR, 12  51 ) ML, IRDl , IRD2 ,  (HMODEL (1,7), 1=1 , 5 ) 
IF(IVSA.EQ.1)CALL  RDNSM 

1251  FORMAT(*0  CARD  2C  ******, 315 , 18A4) 

ENDIF 

ENDIF 

M=7 

CALL  AERNSM ( JPRT ,  GNDALT ) 

IF(ICLD  .LT.  20)  GO  TO  260 
C 

C  SET  UP  CIRRUS  MODEL 

C 

IF (CTHIK. NE.O)  IFLGT=1 
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IF ( CALT . NE . 0 )  IFLGA=1 
IF(ISEED.EQ.O)  IFLGT=2 
IF(ISEED.EQ.O)  I FLG 2 

CALL  CIRRUS  ( CTHIK ,  CALT ,  ISEED ,  CPROB ,  CEXT) 

WRITE (IPR, 1220) 

1220  FORMAT (15X, 'CIRRUS  ATTENUATION  INCLUDED  (N  O  A  A  CIRRUS)  *) 
IF(IFLGT.EQ.O)  WRITE (IPR, 1221)  CTHIK 

1221  FORMAT (15X, 'CIRRUS  ATTENUTION  STATISTICALLY  DETERMENED  TO  BE', 

1  F10.3,'KM*) 

IF(IFLGT.EQ. 1)  WRITE ( IPR, 1222)  CTHIK 

1222  FORMAT  (15X, 'CIRRUS  THICKNESS  USER  DETERMINED  TO  BE  '  ,  FIO .  3  ,  'KM ' ) 
IF(IFLGT.EQ.2)  WRITE (IPR, 1223 )  CTHIK 

1223  FORMAT (15X, 'CIRRUS  THICKNESS  DEFAULTED  TO  MEAN  VALUE  OF  *, 

1  FIO. 3, 'KM') 

I F ( I FLGA . EQ . 0 )  WRITE ( IPR , 1 2  2  4 ) CALT 

1224  FORMAT (15X,  'CIRRUS  BASE  ALTITUDE  STATISCALLY  DETERMINED  TO  BE', 
1  FIO. 3,  *  KM' ) 

IF(IFLGA.EQ.l)  WRITE (IPR, 1225)  CALT 

1225  FORMAT (15X, 'CIRRUS  BASE  ALTITUDE  USER  DETERMINED  TO  BE', 

1  FIO. 3,  *  KM*) 

IF(IFLGA.EQ.2)  WRITE (IPR, 1226)  CALT 

1226  FORMAT  (15X,  'CIRRUS  BASE  ALTITUDE  DEFAULTED  TO  MEAN  VALUE  OF', 

1  F10.3,*KM') 

WRITE ( IPR ,1227) CPROB 

1227  FORMAT ( 15X, 'PROBABILTY  OF  CLOUD  OCCURRING  IS*,F7.1,*  PERCENT') 

C 

C  END  OF  CIRRUS  MODEL  SET  UP 

C 

260  CONTINUE 
C 

c 

C*****CARD  2E 
C 

IF((IHAZE.EQ. 7) .OR. (ICLD.EQ.il))  THEN 
C 

C*****  CARD  2E  USER  SUPPLIED  AEROSOL  EXTINCTION,  ABSORPTION,  AND 

C  ASYMMETRY 

CALL  RDEXA 
C 

ENDIF 

300  CONTINUE 

IPARM  =-99 
IPH  =-99 
IDAY  =-99 
ISOURC=-99 
C 

PARMl  =-99. 

PARM2  =-99. 

PARM3  =-99. 

PARM4  =-99. 

GMT  =-99, 

PSIPO  =  0. 

ANGLEM=-99 . 

G  =-99. 

C 

C*****CARD  3  GEOMETERY  PARAMETERS 
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o  o  o 


driv3210 

IF  ((SEA)  -AND.  (Msea  .EQ.  0))  THEN 
the  first  sky  card  is  next. 

Set  emissivity  to  zero  (TBOUND  is  already  zero)  and  ITYPE  to  3 
for  calculations  coining  up  with  cards  in  TAPES. SEA; 

ITYPE  =  3 
SALE  =0.0 
END  IF 

***  Set  mode  to  sun  irradiance  (if  a  sun  card  will  be  next)  ********* 

IF  ((Sea)  .AND.  (lEMSCTold  .EQ.  2)  .AND.  (Msea  .EQ.  3))  THEN 
lEMSCT  =  3 
IMULT  =  0 
LFIRST  =  .TRUE. 

LSAME  =  .FALSE. 

END  IF 

IF  (lEMSCT  .EQ.  3)  GO  TO  315  driv3220 

*****  Read  introduced  geometry  cards  from  file  TAPE5.SEA  *************** 

IF  (SEA)  THEN 

READ ( IRDS ,1312)H1,H2, ANGLE , RANGE , BETA , RO , LEN , Ps i , SeaSwitch 
ELSE 

READ  ( IRD , 13 12 ) HI , H2 , ANGLE , RANGE , BETA , RO , LEN , Psi , SeaSwitch 

END  IF 

*****  and  remove  the  boundary  (in  sea  AND  sky)  for  a  sea  calculation  *** 

IF  (SeaSwitch)  TBOUND  =0.1 

1312  FORMAT(6F10.3,I5,F10.3,L5) 

WRITE  ( IPR ,  13 13 )  HI ,  H2  ,  ANGLE ,  RANGE ,  BETA ,  RO ,  LEN ,  Psi ,  SeaSwitch 

1313  FORMAT(*0  CARD  3  *****  *  ,  6F10 . 3 , 15  ,  FIO .  3  ,  L5) 

GO  TO  320 

C 

C*****CARD  3  FOR  DIRECTLY  TRANSMITTED  SOLAR  RADIANCE  (lEMSCT  =  3) 

C 

315  CONTINUE 

*********  Read  introduced  sun  card  from  file  TAPE5.SEA  ************** 

IF  (Sea)  THEN 

READ (IRDS, 1316)  H1,H2 , ANGLE, IDAY,RO, ISOURC,ANGLEM 

ELSE 

READ(IRD,  1316)  Hl,H2,ANGLE,IDAY,RO, ISOURC,ANGLEM 

END  IF 

********************************************************************* 

1316  FORMAT(3F10.3,I5,5X,F10.3,I5,F10.3) 

WRITE  ( IPR ,1317)  HI ,  H2  ,  ANGLE ,  IDAY ,  RO ,  ISOURC ,  ANGLEM 

1317  FORMAT(*0  CARD  3  *****  *  ,  3F10 . 3 , 15, 5X,  FIO,  3 , 15  ,F10. 3) 

ITYPE  =  3 
RANGE  =0.0 
BETA  =0.0 
LEN  =  0 

C*****RO  IS  THE  RADIUS  OF  THE  EARTH 
320  RE=6371.23 

C  *********  ERRATA  JULY  25 

IF(H1.  LT.  ZM(1)  )  THEN 

WRITE (IPR, 905)  H1,ZM(1) 

905  FORMAT (*  HI  LESS  THAN  FIRST  ALT  RESET  *,/ 

X  *  HI  WAS  *,F10.2,*  1ST  ALT  =  »,F10.2) 

HI  =  ZM(1) 

ENDIF 

C  *********  end  ERRATA 

HIS  =  HI 
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H2S  =  H2 
ANGLES  =  ANGLE 
RANGS  =  RANGE 
BETAS  =  BETA 
ITYPES  =ITYPE 
LENS  =  LEN 
IF  (MODEL- EQ.O)  RO  =  RE 
IF  (MODEL. EQ.l)  RE=6378.39 
IF  (MODEL. EQ-4)  RE=6356.91 
IF  (MODEL. EQ-5)  RE=6356.91 
IF  (RO.GT.0.0)  RE=RO 
C 

IF  (IEMSCT.NE.2)  GO  TO  330 
C 

C*****CARD  3A1 
C 

IF  (SEA)  THEN 

READ ( IRDS ,1320)  IPARM , IPH , IDAY , I SOURC 

ELSE 

READ ( IRD , 13  2  0 )  IPARM , IPH , IDAY , ISOURC 

END  IF 

1320  FORMAT(4I5) 

WRITE (IPR, 1321)  IPARM, IPH, IDAY , ISOURC 

1321  FORMAT(*0  CARD  3A1***** » , 415) 

C 

C*****CARD  3A2 
C 
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IF  (SEA)  THEN 

READ (IRDS, 1322)  PARMl,  PARM2,  PARM3,  PARM4 , 
+  GMT , PS IPO , ANGLEM , G 

ELSE 

READ (IRD, 1322)  PARMl,  PARM2 ,  PARM3 ,  PARM4 , 
+  GMT , P S IPO , ANGLEM , G 

END  IF 


1322  FORMAT (8F10. 3)  driv3740 

WRITE  ( IPR,  1323)  PARMl ,  PARM2  ,  PARM3 ,  PARM4  ,  GMT ,  PSIPO ,  ANGLEM ,  G 

1323  FORMAT(*0  CARD  3A2***** * , 8F10 . 3 )  driv3760 

C  driv3770 

CSSISSISSISSISSI  CHANGES  BEGIN.  driv3780 

c  driv3790 

REWIND (ISCRCH)  driv3800 

C  driv3810 

IF  (LFIRST  .AND.  IMULT  .EQ.  1)  THEN  driv3820 

C  driv3830 

C  SAVE  SOLAR  PARAMETERS  FOR  COMPARING  LATER,  driv3840 


C  NOTE  THAT  LFIRST  IS  TRUE  AND  IMULT  (MULTIPLE  SOLAR  SCATTERING)  driv3850 
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LFIRST  =  .FALSE. 

CALL  S VSOLA  ( IPARM ,  IPH ,  I  DAY ,  I  SOURC ,  PARMl ,  P ARM2  ,  P ARM3  ,  PARM4  , 
$  GMT, PSIPO, ANGLEM, 

$  ISAVEl ,  ISAVE2  ,  ISAVE3  ,  ISAVE4 ,  SAVEl ,  SAVE2  ,  SAVE3  ,  SAVE4  , 

$  SAVE5  ,  SAVE6  ,  SAVE7 ) 

LSAME  =  .FALSE. 

C 

ELSEIF  (IMULT  .EQ.  1  .AND.  IRPT  .EQ.  3)  THEN 
C 

C  NOW  COMPARE  SOLAR  PARAMETERS;  LSAME  IS  TRUE  IF  THEY  MATCH. 

CALL  COMPAR ( IPARM,  IPH,  IDAY ,  ISOURC ,  PARMl ,  PARM2  ,  PARM3  ,  PARM4 , 
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$  GMT,PSIPO,ANGLEM, 

$  ISAVEl ,  ISAVE2  ,  ISAVE3  ,  ISAVE4  ,  SAVEl ,  SAVE2  ,  SAVE3  ,  SAVE4  , 

$  SAVES ,  SAVE6 ,  SAVE?  ,  LSAME ) 

CALL  SVSOLA  ( IPARM ,  IPH ,  IDAY ,  ISOURC ,  PARMl ,  PARM2  ,  PARM3  ,  PARM4 , 

$  GMT,PSIPO,ANGLEM, 

$  ISAVEl ,  ISAVE2  ,  ISAVE3  ,  ISAVE4  ,  SAVEl ,  SAVE2  ,  SAVES  ,  SAVE4  , 

$  SAVES , SAVE6 , SAVE? ) 

ELSE 

C 

C  GET  READY  FOR  ANOTHER  POSSIBLE  FORTHCOMING  SERIES  OF  MULTIPLE 

C  SOLAR  SCATTERING  RUNS. 

LFIRST  =  .TRUE. 

LSAME  =  .FALSE. 

ENDIF 

c 

CSSISSISSISSISSI  CHANGES  END 
C 

IF{IPH.  EQ  .  0)  THEN 

IF(G.  GE.  1.0)  G  =  .9999 

IF(G.  LE.  -1.0)  G  =  -.9999 

ENDIF 

IF  (IPH.NE.l)  GO  TO  330 
C 

C*****CARD  3B1  USER  DEFINED  PHASE  FUNCTION 
C 

C*****READ  USER  DEFINED  PHASE  FUNCTION 
C 

READ ( IRD ,1326) NANGLS 
1326  FORMAT(IS) 

WRITE ( IPR , 1 3  2  ? ) NANGLS 
132?  FORMAT(*  CARD  3B1***** * , IS) 

C 

C*****CARD  3B2 
C 

READ  (IRD,  1328)  (ANGF{I)  ,  F (1 , 1)  ,  F  (2 , 1)  ,  F  (3 , 1)  ,  F  (4 , 1)  ,1=1, NANGLS) 

1328  FORMAT (5E10. 3) 

WRITE (IPR, 1329)  (ANGF(I) ,F(1,I) ,F(2,I) ,F(3,I)  ,F(4,I)  ,1=1, NANGLS) 

1329  FORMAT(*0  CARD  3B2***** » , 5E10. 3) 

C 

330  CONTINUE 
C 

IF  (IRPT  .EQ.  3)  THEN 

IF  (IPARM  .EQ.  1)  CALL  SUBSOL  (PARM3  ,PARM4  ,  GMT,  IDAY) 

GO  TO  555 
END  IF 
C 

C*****CARD  4  WAVENUMBER 
C 

400  CONTINUE 

IF  (.NOT.  SEA)  THEN 

READ (IRD, * (5110) ' ) IVl , IV2 , IDV, IFWHM, IFILTER 
END  IF 

401  WRITE  (IPR,  *  (15H0  CARD  4  *****,  5110)  *)  IVl ,  IV2  ,  IDV,  IFWHM,  IFILTER 

IF (IDV  .  LE.  0)THEN 

PRINT  ERROR  IN  IDV  »,IDV 

IDV  =  1 
ENDIF 
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IF(IFWHM  .  LE.  0)THEN 

PRINT  ERROR  IN  IFWHM  »  ,  IFWHM 

IFWHM  =  2 
ENDIF 

IF  ((IFILTER  .GE.  1)  .AND.  (IFILTER  .LE.  6))  THEN 
C  reset  wavenumbers  to  span  filter  passband: 

W1  =  FLIST(1,  IFILTER) 

W2  =  FLIST(2,  IFILTER) 

IVl  =  INT(1E4/W2)  -  IDV 
IV2  =  INT(1E4/W1)  +  IDV 

ELSE 

C  filter  data  are  absent.  Reset  to  no  filter  at  all: 

IFILTER  =  0 

END  IF 

IF  (SeaSwitdh)  THEN 

C  Check  number  of  wavenumber  steps.  Reset,  if  necessary,  to 

C  prevent  sea  arrays  from  overflowing  in  "TRANS". 

Nv  =  (IV2  -  IVl) /IDV 

IF  (NV  .GE.  Kv)  IDV  =  (IV2  -  IVl) /Kv  +  1 
END  IF 

WRITE(IP4,  *(1H\,  T20,  22HOUTPUT  FILE  FOR  FILTER,  12, 

+  2H:  ,  15,  3H  TO,  15,  9H  CM“1  IN  ,  12, 

+  12H  CM-1  STEPS.,  /1H\)’)  IFILTER,  IVl,  IV2 ,  IDV 


WRITE (IP4, 

*(1H\,  T65, 

18H  FILTERED 

RADIANCE) 

') 

WRITE (IP4, 

*  (45H\ 

ELEV. 

ANGLE 

RANGE 

TRANS  , 

+ 

T49 

T88 

,  34H  PATH 
,  15H  TOTAL 

SEA 

.  TEMP.)*) 

SKY 

SUN, 

WRITE (IP4, 

* (45H\ 

(mrad) 

(deg) 

(km) 

(--)  , 

+ 

T68 

, 13H  (W  m-2 

sr-1),  TlOO, 

3H(C),  /) 

') 

IF(IHAZE.EQ-3)  THEN 

C  IF(V1.LT. 250.0  .OR.  V2.LT. 250.0)  THEN 

IF(IV1.LT.250)THEN 
IHA2E=4 

WRITE  (IPR,1203) 

ENDIF 

1203  FORMAT  (  *  0**WARNING**  NAVY  MODEL  IS  NOT  USABLE  BELOW  250CM- 

1  /,10X,  *  PROGRAM  WILL  SWITCH  TO  IHAZE=4  LOWTRAN  5  MARITIME* 

END  IF 

IF  (IRPT.EQ.4)  GO  TO  550 
CC  IF  (IRPT.EQ.-4)  GO  TO  560 
500  CONTINUE 

IF  (IRPT.EQ.3)  GO  TO  555 

WRITE (IPR, 1410)  (HTRRAD(I1,IEMSCT+1) ,11=1,6) 

1410  FORMAT(*0  PROGRAM  WILL  COMPUTE  *,6A4) 

IF{ISOURC  .EQ.  1)  WRITE (IPR, 1204) 

1204  FORMAT  (*  LUNAR  SOURCE  ONLY  ') 

IF  (IMULT  .EQ.  1)  THEN 

IF(IEMSCT.EQ.  0  .OR.  IEMSCT.EQ.3  )  THEN 
WRITE (IPR, 1411) 

1411  FORMAT(»0  MULTIPLE  SCATTERING  HAS  BEEN  TURNED  OFF  *) 
WRITE  (IP6, 

+  *(/,  39HMULTIPLE  SCATTERING  HAS  BEEN  TURNED  OFF)*) 

IMULT=0 

ELSE 

WRITE(IPR,1412) 

END  IF 
END  IF 
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1412  FORMAT{*0  CALCULATIONS  WILL  BE  DONE  USING  MULTIPLE  SCATTERING  *) 
MDEL=MODEL 
IF (MDEL . EQ . 0 ) MDEL=8 
MM1=MDEL 
MM2=MDEL 
MM3=MDEL 

IF(M1.NE.0)MM1=M1 
IF(M2.NE.0)MM2=M2 
IF(M3,NE.0)MM3=M3 
IF (MODEL. EQ.O)  GO  TO  510 

WRITE  (IPR;  1500)  MMl,  {HM0DEL(I1,MM1)  ,11=1,5)  ,MM2,  (HMODEL(I2  ,MM2 )  , 
1  12=1,5) ,MM3,  (HMODEL(I3,MM3) ,13=1,5) 

1500  FORMAT ( ’ 0  ATMOSPHERIC  MODEL  * , /  , 

1  lOX, ‘TEMPERATURE  =  * , 14 , 5X, 5A4  ,  /  , 

1  lOX, ‘WATER  VAPOR  =  * , 14 , 5X, 5A4 , /  , 

1  lOX, ‘OZONE  =  »,I4,5X,5A4) 

WRITE (IPR, 1501)  M4,M5,M6,MDEF 

1501  FORMAT(20X,’  M4  =  ‘,15,  ‘  M5  =  ',15,*  M6  =  *,I5,*  MDEF  =  *  ,15) 

C 

510  IF(JPRT.EQ.O)  GO  TO  520 
IF ( ISEASN . EQ . 0 ) I SEASN=1 
IF(IVULCN.LE.O)  IVULCN=1 
IHVUL=IVULCN+10 
IF(  IVULCN  .EQ.  6)  IHVUL  =  11 
IF(  IVULCN  .EQ.  7)  IHVUL  =  11 
IF(  IVULCN  .EQ.  8)  IHVUL  =  13 
IHMET=1 

IF ( IVULCN - GT . 1 ) IHMET=2 
IF(IHAZE.EQ.O)  GO  TO  520 

WRITE  (IPR,  1510)  (HHAZE(I,IHAZE)  ,1=1,5)  ,VIS,  (HHAZE(I2,6)  ,12=1,5)  , 

1  (HHAZE(II,6)  ,11=1,5)  ,  (HSEASN(IA,  ISEASN)  ,IA=1, 5)  , 

2  (HHAZE(I3, IHVUL)  ,13=1,5)  , 

3  (HVULCN( IB,  IVULCN)  ,IB=1,5)  ,  (HSEASN (IC,  ISEASN)  ,IC=1,5)  , 

4  (HHAZE(I4, 16)  ,14=1,5)  ,  (HMET(I5,IHMET)  ,15=1,5) 

1510  FORMAT(»0  AEROSOL  MODEL  *,/,  lOX, ‘REGIME*  , 

A  T35,  ‘AEROSOL  TYPE»,T60,  ‘PROFILE* ,T8 5,  ‘SEASON*,/,/, 

B  1 OX, ‘BOUNDARY  LAYER  (0-2  KM) * , T35 , 5A4 ,T60 , F5 . 1, 

C  ‘  KM  VIS  AT  SEA  LEVEL  *,/,  lOX,  *  TROPOSPHERE  (2-lOKM)  ‘  ,  T35 , 

D  5A4,T60,5A4,T85,5A4,/,10X,  ‘STRATOSPHERE  (10-30KM)  *  , 

E  T35,5A4,T60,5A4,T85,5A4,/,10X, 'UPPER  ATMOS  (30-100KM) ’ , 

F  T35,5A4,T60,5A4) 

520  CONTINUE 

IF(ITYPE.EQ.l)  THEN 

WRITE(IPR,1515)  HI, RANGE 
WRITE(IP6,1515)  HI, RANGE 
END  IF 

1515  FORMAT ( / ,  *  HORI ZONTAL  PATH  * , / /  , 

1  8X, ‘ALTITUDE  =  ‘,F10.3,'  KM*,/, 

2  8X, ‘RANGE  =  ‘,F10.3,‘  KM‘,/) 

IF(ITYPE.EQ.2)  THEN 

WRITE(IPR, 1516)  H1,H2, ANGLE, RANGE, BETA, LEN 
WRITE  ( IP6 ,1516)  HI ,  H2  ,  ANGLE ,  RANGE ,  BETA ,  LEN 
END  IF 

1516  FORMAT (/,*  SLANT  PATH,  HI  TO  H2*,//, 

1  10X,*H1  =  ‘,F10.3,*  KM‘,/,10X,  •H2  =  *,F10.3,’  KM*,/, 

2  lOX,  ‘ANGLE  =  ',F10.3,‘  DEG*  ,/, lOX, ‘RANGE  =  ‘,F10.3,‘  KM‘,/, 

3  10X,‘BETA  =  *,F10.3,‘  DEG  *  ,  /  ,  lOX,  *  LEN  =  *,I6,/) 
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IF(ITYPE.EQ.3)  THEN 

WRITE (IPR, 1517)  HI, H2, ANGLE 
WRITE (IP6, 1517)  HI, H2, ANGLE 
END  IF 


1517 

FORMAT (/ 

,  *  SLANT 

PATH  TO  SPACE',//, 

1 

lOX,  ' 

HI 

=  • ,F10.3, ' 

KM’,/, 

2 

lOX,  ' 

HMIN 

=  ' ,F10.3, ' 

KM',/, 

3 

lOX,  ' 

ANGLE 

=  ’,F10.3,' 

DEG',/) 

IF  (IEMSCT.NE.2)  GO  TO  550 
C 

C*****INTREPRET  SOLAR  SCATTERING  PARAMETERS 

C 

C 

IF  (IPARM.EQ.l)  CALL  SUBSOL  (PARM3  ,  PARM4  ,  GMT,  IDAY) 


1530 


1532 


WRITE  (IPR,1530) 

FORMAT(’0  SINGLE  SCATTERING  CONTROL  PARAMETERS  SUMMARY  */) 
IF(IPARM.NE.2)  WRITE  (IPR,1532)  PARMl , PARM2  , PARM3  , PARM4 , GMT , PSIPO 
1,IDAY 

FORMAT (1 OX, 'OBSERVER  LATITUDE  =' ,T35,F10.2,  '  DEG  NORTH  OF  EQUATOR' 

1  ,/,10X, 'OBSERVER  LONGITUDE= * ,T35 , FIO .  2 ,  *  DEG  WEST  OF  GREENWICH', 

2  /,10X,  'SUBSOLAR  LATITUDE  = » , T35 , FIO .  2  ,  '  NORTH  OF  EQUATOR*,/, 

3  lOX,  'SUBSOLAR  LONGITUDE  = * ,T35 , FIO .  2 ,  '  WEST  OF  GREENWICH*,/, 

4  10X,'TIME  (<0  IS  UNDEF)  =  ' ,T35,F10.3,  *  GREENWICH  TIME*,/, 

5  1 OX, 'PATH  AZIMUTH  = ' , T35 , FIO • 3 , *  DEG  EAST  OF  NORTH*,/, 

6  10X,'DAY  OF  YEAR  =',T35,I10) 

IF  (IPARM.EQ.2)  WRITE  (IPR,  1534 )  PARMl , PARM2  , GMT, PSIPO,  IDAY 
FORMAT (1 OX, 'RELATIVE  AZIMUTH  = * , T35, FIO . 3  ,  *  DEG  EAST  OF  NORTH*,/, 

1  lOX, 'SOLAR  ZENITH  = ' ,T35 , FIO . 3 , '  DEG  *,/, 

2  10X,»TIME  (<0  UNDEF)  = » , T35 , FIO . 3 ,  '  GREENWICH  TIME ',/ , 

3  10X,'PATH  AZIMUTH  = * , T35 , FIO . 3 , '  DEG  EAST  OF  NORTH*,/, 

4  10X,'DAY  OF  THE  YEAR  =',T35,I6) 

IF  (ISOURC.EQ.O)  WRITE  (IPR, 1535) 

FORMAT('0  EXTRATERRESTIAL  SOURCE  IS  THE  SUN') 

IF  (ISOURC.EQ.l)  WRITE  (IPR, 1536)  ANGLEM 

FORMAT('0  EXTRATERRESTIAL  SOURCE  IS  THE  MOON,  MOON  PHASE  ANGLE  =*, 
1  FIO- 2, »  DEG*) 

IF  (IPH.EQ.O)  WRITE  (IPR, 1538)  G 
F0RMAT('0  H-G  PHASE  FUNCTION  ,G=*,F10-3) 

IF  (IPH.EQ.l)  WRITE  (IPR, 1540) 

FORMAT(*0  USER  SUPPLIED  PHASE  FUNCTION*) 

IF  (IPH.EQ-2)  WRITE  (IPR, 1542) 

FORMAT(*0  PHASE  FUNCTION  FROM  MIE  DATA  BASE') 

CONTINUE 

VI  =FLOAT(INT{V1/5.0+0.1) )*5-0 
V2  ==FLOAT(INT(V2/5.0+0-1))*5.0 
TO  AVOID  THE  DIFFICULTY  FOR  V1=0 
ALAM1=  99999.98 
IF (VI . GT . 0 . ) ALAM1=10000 . /VI 
ALAM2=10000./V2 
IF(DV.LT.5-)DV=5. 

DV=FLOAT(INT(DV/5+0.1) ) *5.0 
WRITE  (IPR, 1555)  VI ,  ALAMl ,  V2 ,  ALAM2  ,  DV 
C1555  FORMAT(*0  FREQUENCY  RANGE  */,10X,*  VI  =  ',F12.1,'  CM-1  (*, 

C  1  FIO. 2,*  MICROMETERS) * ,/,l OX, '  V2  =  *,F12.1,'  CM-1  (*,F10.2, 

C  2  '  MICROMETERS) ' ,/lOX, *  DV=  *,F12.1,*  CM-1*) 

IF ( . NOT . MODTRN) THEN 
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C 
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IVl=5*(IVl/5) 

IV2=5*(  (IV2-f4)  /5) 

IDV=5+5*((IDV-5)/5) 

ENDIF 

IF ( IV2 . LT . IVl+IDV) THEN 

WRITE  (IPR,  *  (/41H  IV2  WAS  LESS  THAN  IVl  4*  IDV  AND  HAS  BEEN, 
1  6H  RESET,/)') 

IV2=IV1+IDV 

ENDIF 


CRZ 

IF(MODTRN)THEN 

CRZ 

IV1SAV=IV1 

CRZ 

IV2SAV=IV2 

CRZ 

IDVSAV=IDV 

CRZ 

ENDIF 

I F ( IVl . NE . 0 ) ALAM1=10  0  00,/ IVl 

ALAM2=10000,/IV2 

IF ( IFWHM . LT . 1 ) IFWHM=1 

IF ( IFWHM . GT . 50 ) IFWHM=5  0 

WRITE (IPR, ' (17H0  FREQUENCY  RANGE, /lOX, 8H  IVl  =,I10,8H  CM-1  ( 

1  F10.2,13H  MICROMETERS) ,/10X,8H  IV2  =,110, 8H  CM-1  (,F10.2, 

2  13H  MICROMETERS) ,/10X,8H  IDV  =,110, 5H  CM-1,/10X,8H  IFWHM  =, 

3  110, 5H  CM-1) ») IVl, ALAMl, I V2,ALAM2, IDV, IFWHM 

WRITE ( IP6, ' (15HFREQUENCY  RANGE, //lOX, 9HIV1  =,I11,8H  CM-1  ( 

1  F7.2,13H  MICROMETERS) ,/10X,9HIV2  =,111, 8H  CM-1  ( , F7 • 2 , 

2  13H  MICROMETERS) ,/10X,9HIDV  =,I11,5H  CM-1, /lOX, 9HIFWHM 

3  I11,5H  CM-1,/10X,9HIFILTER  =,111) ») 

4  IVl ,  ALAMl ,  IV2  ,  ALAM2  ,  IDV,  IFWHM,  I  FILTER 
C 

C*****lOAD  ATMOSPHERIC  PROFILE  INTO  /MODEL/ 

C 

CALL  STDMDL 
C 

C  DEFINE  COUNTER  ITEST  TO  PREVENT  ZENITH  ANGLE  QTHETA  AND  LAYER 
C  PATH  LENGTH  PL  FROM  BEING  CHANGED  DURING  SOLAR  CALCULATIONS 
555  DO  15  1=1,102 

DO  15  J=1,KMAX 

WPATH(I, J)=0.0 
15  WPATHS(I, J)=0,0 
C 

ITEST=0 

C 

IF  (IMULT  .EQ.  1)  THEN 
H1=ZM(1) 

H2=ZM(ML) 

ITYPE  =  2 
ANGLE  =  0. 

BETA  =  0. 

RANGE  =0. 

ISSGS  =  ISSGEO 
ISSGEO  =  0 

C  CALL  GEO  (IERROR,BENDNG,MAXGEO) 

MSOFF=68 

CALL  GEO  (IERROR,BENDNG,MAXGEO,MSOFF) 

W15SV  =  W(15) 

C 

C  W15SV  IS  THE  REL  HUM  FROM  0  TO  SPACE 

C  THIS  REL  HUM  MAY  BE  DIFFERENT  THAN  THE  PATH  REL  HUM 


driv5900 

driv5910 

driv5920 

driv5930 

driv5940 

driv5950 

driv5960 

driv5970 

driv5980 

driv5990 

driveooo 

driveolo 

driv6020 

driv6030 

driv6040 

driv6050 

driv6060 

driv6070 

driv6080 

driv6090 

driveiOO 

drivSllO 


driv6120 
driv6130 
driv6140 
driv6150 
driv6160 
drive 170 
driv6180 
driv6190 
driv6200 
driv6210 
driv6220 
driv6230 
driv6240 
driv6250 
driv6260 
driv6270 
driv6280 
driv6290 
driv6300 
drives 10 
driv6320 
drives 30 
o  driv6340 
driv6350 
driveseo 

driv6370 

driv6380 

driv6390 

driv6400 

driv6410 


B-19 


*  n  n  *  no  noon  n  nooo  no 


WHEN  REL  HUM  ARE  DIFFERENT  THE  ANSWER  CAN  CHANGE 

ISSGEO  =  ISSGS 
IMSMX=IKMAX 
DO  35  N=1,IMSMX 
PLST(N)=PL(N) 

DO  35  K==1,KMAX 
5  WPMS ( N , K ) =WPATH ( N , K ) 

35  PLST(N)=PL(N) 

IF ( lEMSCT • EQ . 2 )  THEN 

CALL  SSGEO  ( TERROR ,  IPH ,  IPARM ,  PARMl ,  PARM2 , 

1  PARM3 , PARM4 , P S I PO , G , MAXGEO ) 

1  PARM3  ,  PARM4  ,  PSIPO ,  G ,  MAXGEO ,  MSOFF ) 

DO  30  N=1,IKMAX 

CSENSV{N)  =  ABS(CSZEN(N) ) 

IF(CSENSV(N)  .LT.  0.0174)  CSENSV(N)  =  0.0174 
30  CONTINUE 

DO  45  N=1,ML 

DO  45  K=1,KMAX 

WPMS  S ( N , K ) =WP ATH S ( N , K ) 

45  CONTINUE 

ENDIF 
ENDIF 

HI  =  HIS 

H2  =  H2S 

ANGLE  =  ANGLES 
RANGE  =  RANGS 
BETA  =  BETAS 
ITYPE  =  ITYPES 
LEN  =  LENS 

****traCE  PATH  THROUGH  THE  ATMOSPHERE  AND  CALCULATE  ABSORBER  AMOUNTS 

ISSGEO=0 

MSOFF=0 

****  Save  original  value  of  SEA  (false  if  earth  not  yet  hit)  ****** 
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SEAold  =  SEA 

CALL  GEO(IERROR,BENDNG, MAXGEO, MSOFF)  driv6790 

***  and  set  HIT  true  if  the  earth  has  been  hit  within  FNDHMN:  ***** 

HIT  =  ((.NOT.  SEAold)  .AND.  SEA) 

IF  (HIT)  THEN 

WDT'P'P  /TPA  *  (  /  fiHC’pa  AT  F7  7 

+  31H  K  REPLACES  BLACK  BODY ' BOUNDARY, //, lOX, 9HUPWIND  =,F11.3, 

+  26H  DEG  EAST  OF  LINE  OF  SIGHT)  *)  TBOUNDold,  Psi 

C  calculate  geometry  from  point  of  view  of  the  footprint 

IF  {lEMSCTold  .EQ.  1)  Pr  =  Psi*d2r  +  pi 
IF  ((IPARM  .EQ.  0)  .OR.  (IPARM  .EQ.  1))  THEN 
ThetaO  =  PARMl 
PhiO  ==  PARM2 
Thetas  =  PARM3 
Phis  =  PARM4 

CALL  Foot (ThetaO, PhiO , ThetaS , PhiS , PsiPO , Beta , Ps i) 

ELSE  IF  (IPARM  .EQ.  2)  THEN 
PsiO  =  PARMl 
DelO  =  PARM2 
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o  o  o 


CALL  SunFoot (PsiO , DelO , PsiPO, Beta , Psi ) 

END  IF 

C  and  issue  new  sky  (and  sun)  cards  in  file  ’TAPE5.SEA* 

CALL  Card 
END  IF 

IF  ((SeaSwitch)  .AND.  (.NOT.  Sea))  THEN 

WRITE  (IP6,  *(/,  13HTBOUND  SET  TO,  F7.2, 

+  17H  K  FOR  MARINE  SKY) »)  TBOUND 

END  IF 
C 

CALL  AERTMP 

IF(IMULT.  NE.  1)  W15SV  =  W(15) 

SAVE  TEMPERATURE  AND  PATH  INFO  FOR  LATER  USE 

IF(IMULT  .EQ,  1)  THEN 
DO  25  N=1,IKMAX 
25  QTHETS(N)  =  QTHETA(N) 

ENDIF 
C 

IF(IERROR.GT.O)  GO  TO  630 

IF(IEMSCT.EQ.3  .AND.  TERROR. EQ.  -5)  GO  TO  557 
GO  TO  558 

557  CONTINUE 
WRITE (IPR, 1557) 

1557  FORMAT(*0  DIRECT  PATH  TO  SUN  INTERSECTS  THE  EARTH:  SKIP  TO 
1  *NEXT  CASE*) 

GO  TO  630 

558  CONTINUE 
C 

IF  ( lEMSCT .  EQ .  2 )  CALL  SSGEO  ( TERROR ,  IPH ,  IPARM ,  PARMl ,  PARM2 ,  PARM3  , 

C  1  PARM4,PSIPO,G,MAXGEO) 

1  PARM4 , PSIPO , G , MAXGEO , MSOFF ) 

W(15)  =  W15SV 
C 

C  W15SV  IS  THE  REL  HUM  (FOR  MULT  SCAT  THIS  MAY  BE  DIFFERENT 

C  FROM  PATH  REL  HUM) 

C 

C  THE  SECOND  CALL  TO  SSGEO  IS  TO  GET  THE  CORRECT  ANGLES  FOR 

C  PHASE  FUNCTIONS 

C 

C  SAVE  SOLAR  PATH  INFORMATION 

C 

IF ( TERROR. GT.O)  GO  TO  630 
C 

IF(IMULT,EQ.l)  THEN 

DO  60  IK  =  1,IMSMX 
PL(IK)=PLST(IK) 

IF ( lEMSCT . EQ . 2 )  CSZEN ( IK) =CSENSV ( IK) 

60  CONTINUE 

DO  70  IK  =  1,IKMAX 

70  QTHETA(IK)  =  QTHETS(IK) 

ENDIF 

C 

C*****L0AD  AEROSOL  EXTINCTION,  ABSORPTION,  AND  ASYMMETRY  COEFFICIENTS 
C 

CALL  EXABIN 
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C  driv7260 

C*****WRITE  HEADER  DATA  TO  TAPE  7  driv7270 

C  driv7280 

C560  WRITE(IPU,1110)MODEL,ITYPE,IEMSCT,IMULT,M1,M2,M3,  driv7290 

C  1  M4,M5,M6,MDEF,  IM,NOPRT,TBOUND,SALB  driv7300 

560  WRITEdPU,  '  (L1,I4,12I5,F8.3,F7.2) ’)M0DTRN,M0DEL  driv7310 

1  ,ITYPE,IEMSCT,IMULT,M1,M2,M3,M4,M5,M6,MDEF,IM,N0PRT,TB0UND,SALB  driv7320 
C  WRITE(IPR1,1110)MODEL,ITYPE,IEMSCT,IMULT,M1,M2,M3,  driv7330 

C  1  M4,M5,M6,MDEF,  IM,N0PRT,TB0UND,SALB  dxiv7340 

WRITE (IPRl, ’ (L1,I4,12I5,F8.3,F7.2) ' )M0DTRN,  MODEL  driv7350 

1  ,  ITYPE ,  IEMSCT  ,  IMULT ,  Ml ,  M2 ,  M3 ,  M4 ,  M5 ,  M6 ,  MDEF,  IM,  NOPRT ,  TBOUND,  SALB  dr iv7  360 


WRITE  (IPU,  1200)  IHAZE,  ISEASN,  IVULCN,  ICSTL,  ICLD,  IVSA,  VIS,WSS,  WHH, 

1  RAINRT,GNDALT 

WRITE  (IPRl,  1200)  IHAZE,  ISEASN,  IVULCN,  ICSTL,  ICLD,  IVSA/  VIS,  WSS,  WHH, 
1  RAINRT,GNDALT 

WRITE (IPU, 1210)  CTHIK, CA1T,CEXT,  ISEED 
WRITE (IPRl, 1210)  CTHIK, CALT, CEXT,  ISEED 
WRITE  ( IPU,  1230)  ZCVSA,  ZTVSA,  ZINVSA 
WRITE  ( I  PRl ,  12  3  0 )  ZCVSA,  ZTVSA,  Z INVSA 
WRITE (IPU, 1255)  ML,  (HMODEL (I, 7) ,  1=1, 5) 

WRITE (IPRl, 1255)  ML,  (HMODEL(I, 7) ,  1=1,  5) 

1255  FORMAT (  I5,18A4) 

IF (MODEL. NE.O) WRITE  (IPU, 1312)  HI, H2, ANGLE, RANGE, BETA, RO, LEN 
IF (MODEL. NE.O) WRITE  (IPRl, 1312)  HI, H2, ANGLE, RANGE, BETA, RO, LEN 
HMDLZ(8)  =  RANGE 
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IF  (MODEL. EQ. 0)  WRITE (IPU, 1560)  (HMDLZ (K) , K-1, 8) 
IF(MODEL.EQ.O)  WRITE(IPR1, 1560) (HMDLZ (K) ,K=1, 8) 

1560  FORMAT(3F10.3, 5E10.3) 

WRITE (IPU, 1320)  IPARM,  IPH, IDAY,  ISOURC 
WRITE (IPRl, 1320)  IPARM,  IPH, IDAY, ISOURC 

WRITE  { I PU ,  1 3 2  2 )  PARMl ,  PARM2 ,  PARM3 ,  PARM4 ,  GMT ,  PS  I PO ,  ANGLEM,  G 
WRITE  (IPRl,  1322)  PARMl,  PARM2,  PARM3,  PABM4,GMT,  PSIPO,  ANGLEM,  G 
3  WRITE (IPU, 1400)  V1,V2,DV 

3  WRITE (IPRl, 1400) VI,  V2,  DV 

WRITE (IPU,  '  (5110) ’ ) IVl, IV2, IDV, IFWHM, IFILTER 
WRITE (IPRl, '  (5110) ' ) IVl, IV2, IDV, IFWHM, IFILTER 


CRZ  IRAIN=0 

CRZ  IF(RAINRT.GT.O)  IRAIN=1 

ccc 

CCC  CALCULATE  EQUIVALENT  LIQUID  WATER  CONSTANTS 
CCC 

CALL  EQULWC 
IF  (SEA)  THEN 

Msea  =  Msea  +1 

Done  =  (((lEMSCTold  ,EQ.  1)  .AND.  (Msea  .EQ.  3)) 
+  .OR.  ((lEMSCTold  .EQ.  2)  .AND.  (Msea  .EQ.  4))) 

IF  (Done)  THEN 

READdRD,  1600)  IRPT 

ELSE 

READ (IRDS, 1600)  IRPT 

END  IF 

ELSE 

READdRD,  1600)  IRPT 

END  IF 

1600  FORMAT (15) 

WRITEdPU,  1600)  IRPT 
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C  TRANS  RETURNS  IRPT  =  -4  IF  THE  SPECTRAL  RANGE  EXTENDS  BEYOND  THE 

C  BAND  MODEL  TAPE.  IN  THIS  CASE,  A  LOWTRAN  7  CALCULATION  IS 

C  PERFORMED  FOR  THE  SHORT  WAVELENGTHS  AND  THEN  THE  ORIGINAL  INPUT 

C  IS  RESTORED.  (NOTE:  THIS  FEATURE  WAS  COMMENTED  OUT.) 

C 

*************  Reset  the  parameters  to  their  original  values  ******** 
*************  from  before  TAPES. SEA  was  introduced,  provided  ******** 
*************  all  cards  from  TAPES. SEA  have  been  read.  ************** 
C 


IF  (Done)  THEN 


CLOSE  (IRDS,  STATUS= ' KEEP ' ) 

Sea  =  . FALSE . 

SeaOld  =  .FALSE. 

Hit  =  .FALSE. 

Msea  =  -1 

ITYPE  =  ITYPEold 

lEMSCT  =  lEMSCTold 

IMULT  =  IMULTold 

TBOUND  =  TBOUHDold 

SALE  =  SALBold 

END  IF 

C 

*********************************************************************** 

C*****WRITE  END  OF  FILE  ON  TAPE  7 

drivSOSO 

630 

IF (TERROR  .GT.  0)  THEN 

driv8060 

READ (IRD, 1600, END=900)  IRPT 

driv8070 

WRITE (IPU, 1600)  IRPT 

drivSOSO 

WRITE ( IPRl ,1600) IRPT 

driv8090 

ENDIF 

drivSlOO 

WRITE(IPU,1620) 

drivSllO 

WRITE (IPRl, 1620) 

driv8120 

1620 

FORMAT (»  -9999.*) 

driv8130 

C 

driv8140 

WRI TE ( IPR , 1 6  3  0 ) I RPT 

drivSlSO 

1630 

FORMAT(*0  CARD  S  ******,15) 

driv8160 

IF  (IRPT.EQ.O)  GO  TO  900 

drive 170 

IF  (IRPT.EQ.4)  GO  TO  400 

driv8180 

cssi 

IF  (IRPT.GT.l  .AND.  IEMSCT.EQ.3)  THEN 

driv8190 

cssi 

PRINTgl,'/!!  ERROR  IN  INPUT  lEMSCT  EQ  3  IRPT  GT  11 » 

driv8200 

cssi 

STOP 

driv8210 

cssi 

ENDIF 

driv8220 

IF  (IRPT.GT.4)  GO  TO  900 

driv8230 

GO  TO  (100,900,300,400),  IRPT 

driv8240 

900 

Call  Timer (lend) 

C 

to  find  how  long  the  calculation  took: 

1880 

WRITE  (ITM,  1880)  ( FLOAT (lend) -FLOAT (Istart) ) /lOO . 

FORMAT  (‘Elapsed  time  (sec)  for  the  last  run  was  *,  F8.2) 

STOP 

driv8250 

END 

driv8260 
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WRITE (IPRl, 1600) IRPT 
C 

ground= .false . 

if (h2. le.zm(l) ) ground^ . true . 

IF  (Msea  .GT.  -1) 

+  PRINT  ’(35H  Driver:  Calling  TRANS  for  sea  card,  12,  IH.) 
CALL  TRANS ( IPH , I SOURC , IDAY , ANGLEM , ground ) 

C 
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APPENDIX  C 

MODIFIED  SUBROUTINE  “TRANS” 


C-1 


N  INTO 
RADSUM 


DONE? 


CALCULATE 
SEA  RADIANCE. 


RETURN 
FROM  TRANS. 


Figure  C-1 .  Flowchart  for  modified  subroutine  “TRANS.” 
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c 

c 

c 


c 


c 


c 

c 


c 

c 

c 

C 


SUBROUTINE  TRANS  (IPH,  ISOURC,  IDAY,  ANGLEM,  ground) 

CALCULATES  TRANSMITTANCE  AND  RADIANCE  VALUES  BETWEEN  IVl  AND  IV2 
FOR  A  GIVEN  ATMOSPHERIC  SLANT  PATH 
parameter  (nbins=99 ,  iprint=50,maxv=50000) 
real  WGT  ( nbins )  ,  SLIT  (56,  nbins ) 

LOGICAL  IVTEST ,  loopO ,  ground ,  transm , modtrn 

COMMON  RELHUM(34)  ,WHN03(34)  ,ICH(4)  ,  VH ( 17 )  ,  TX (63 )  ,  W (63  )  ,IMSMX, 

1  WPATH(102,63)  ,TBBY(102)  ,PATM(102)  ,NSPEC,KP0INT(12)  , 

2  ABSC(5,47) ,EXTC(5,47) ,ASYM(5,47) ,VX0(47) ,AWCCON(5) 

COMMON  /IFIL/  IRD,IPR,IPU,NPR,IPR1,IP6,IP7,IP8,IP4,IRDS,IP6S, 

+  ITR, Isky, Isun,Ipath 

COMMON / CARD  1  /MODEL ,  ITYPE ,  lEMSCT ,  Ml ,  M2 ,  M3  ,  IM ,  NOPRNT ,  TBOUND ,  SALE , 

1  MODTRN 

COMMON  /CARD2/  IHAZE,  ISEASN,  IVULCN,  ICSTL,  ICIR,  IVSA,VIS,  WSS,  WHH, 

1  RAINRT 

COMMON/ CARD3  /  HI ,  H2  ,  ANGLE ,  RANGE ,  BETA ,  REE ,  LEN 
COMMON/CARD4/IV1,  IV2  ,  IDV,  IFWHM,  IFILTER 
COMMON /  CNSTNS /PIX ,  CA ,  DEG ,  GCAIR ,  BIGNUM ,  BIGEXP 
COMMON / CNTRL/  KMAX ,  M ,  IKMAX ,  NL ,  ML ,  IKLO ,  I SSGEO ,  IMULT 
COMMON/ SOLS /AH1( 68)  ,ARH(68)  , 

1  WPATHS(102,63)  ,PA(68)  ,PRX(68)  ,ATHETA(35)  ,ADBETA(35)  ,LJ(69)  , 

2  JTURN,ANGSUN,CSZEN(68)  ,TBBYS  (102 , 12)  ,PATMS  (102 , 12) 
COMMON/SRAD/TEBl , TEB2SV 

COMMON/MSRD/TLE(34)  ,COSBAR(34)  ,OMEGAO(68)  ,UPF(10,34)  ,DNF(10,34)  , 

1  TAER(34)  ,ASYIK(68)  ,ASYDM(68)  ,STRN(0:34)  ,DMOLS(68)  ,DSTRN(0:68)  , 

2  FDNSRT,FDNTRT,TAUT(34)  ,UMF(34)  ,DMF(34)  ,UMFS(34)  ,DMFS(34) 
COMMON/ICLL/ICALL,  FPHS ,  FALB,  FORBIT 

PARAMETER  (Kr  =  216,  Kv  =  400) 

COMMON  /Filters/  FLIST(5,6), 

+  FILTER1(45),  BBl(Kr),  FILTER2(54),  BB2 (Kr) , 

+  FILTER3(39),  BB3 (Kr) ,  FILTER4(47),  BB4 (Kr) , 

+  FILTERS ( 101) , BBS (Kr) ,  FILTER6(75),  BB6(Kr) 

COMMON/ Sea /  Sea ,  Hit ,  Msea ,  TBOUNDold ,  lEMSCTold 
COMMON /Geometry/  To,Po,Tr,Pr 

COMMON / Constants /pi ,  r2 d ,  d2 r ,  epsilon ,  delta ,  onem ,  onep ,  inf  inity 
DIMENSION  Tau(Kv),  SkyN(3,Kv),  Ho(Kv) ,  Tsky(3)  ,  Rsky(3) 

LOGICAL  Sea,  PathCard,  SkyCard,  LastSky,  SunCard,  Hit 

REAL  Npath,  Nsky,  Nvsky,  Nsea,  Nsun,  Ntotal,  Nbb,  infinity.  No 

TO  =  273.15 

common  /solar/ Isame 

logical  Isame 

Initialize  slit  function  array 
DO  10  I  =  1,56 

DO  10  J  =  1, nbins 
10  SLIT(I,J)  =0. 

Initialize  radiance  minimum  and  maximum  parameters 

RADMIN=bignum 

RADMAX=0 . 

Initialize  ground  emissivity  (one  minus  ground  albedo) 
EMISS=1.“SALB 

Store  the  number  of  path  layers  in  ikmx 
IKMX==IKMAX 


tras  110 
tras  120 
tras  130 
tras  140 
tras  150 
tras  160 
tras  170 
tras  180 
tras  190 


tras  210 
tras  220 


tras  250 
tras  260 

tras  280 
tras  290 
tras  300 
tras  310 
tras  320 
tras  330 


tras  340 
tras  350 
tras  360 
tras  370 
tras  380 
tras  390 
tras  400 
tras  410 
tras  420 
tras  430 
tras  440 
tras  450 
tras  460 
tras  470 
tras  480 
tras  490 
tras  500 
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noon 


c 

c 

c 


c 

c 

c 

c 

c 

c 

cssi 


tras 

510 

Initialize 

integrated  absorption,  radiance,  solar  irradiance  and 

tras 

520 

transmitted 

■  solar  irradiance  sums 

tras 

530 

SUMA=0. 

tras 

540 

RADSUM=0 

tras 

550 

SSOL=0 . 

tras 

560 

STSOL=0 . 

tras 

570 

Path Sum 

= 

0. 

PathCard 

= 

.FALSE. 

SkyCard 

.FALSE. 

LastSky 

.FALSE, 

SunCard 

= 

.FALSE. 

IF  (Sea) 

THEN 

IF  (Msea  .EQ.  0)  THEN 

PathCard  =  .TRUE. 


ELSE  IF  ((Msea  .GE.  1)  .AND.  (Msea  .LE.  3))  THEN 
SkyCard  =  .TRUE. 

ELSE  IF  (Msea  .EQ.  4)  THEN 
SunCard  =  .TRUE. 

END  IF 

IF  (Msea  .EQ.  3)  THEN 
LastSky  =  .TRUE. 

END  IF 
END  IF 

IF  (SkyCard)  Tsky(Msea)  =  ANGLE*d2r 
Istore  =  0 

tras  580 


Initialize  integration  weighting  factor  tras  590 

FACT0R=.5  tras  600 

tras  610 

Initialize  icount,  used  to  determine  when  header  must  be  printed  tras  620 

ICOUNT=iprint  tras  630 

tras  640 


Do  not  perform  a  MODTRAN  calculation 
IF ( IVl . GE. 2  2  655 ) modtrn= . false . 

IF ( IVl . GE . 22  68 1 ) modtrn= .false . 
IF(modtrn)THEN 


if  all  sources  are  continuum  tras  650 

tras  660 
tras  670 
tras  680 
tras  690 


WHEN  THE  band  model  or  line-by-line  option  is  used,  call 
routine  "bmdata”  to  INITIALIZE  PARAMETERS  AND  TO  SET  THE 
FREQUENCY  STEP  SIZE  **IDVX"  TO  THE  BAND  WIDTH  (1  CM-1)  . 
IDV5=5 

CALL  bmdata  ( IVl ,  IFWHM ,  IDVX ,  IKMX ,  MXFREQ ) 

IWIDM1=IFWHM/ IDVX-1 
IV=5* ( (IVl-IWIDMl) /5) 

IF(IV.LT,0)IV=0 

IVX=IV-IDVX 

IV=IV-5 

IVXMAX=IV2+IWIDM1 


tras  700 
tras  710 
tras  720 
tras  730 
tras  740 
tras  750 
tras  760 
tras  770 
tras  780 
tras  790 
tras  800 


ELSE 

IDV5=IDV 

IDVX=IDV5 

IWIDM1=0 

IV=IV1-IDV5 

IVX=IV 

IVXMAX=IV2+IWIDM1 

IF  ( IVXMAX .  GT .  maxv)  IVXMAX^maxv 


tras  810 
tras  820 
tras  830 
tras  840 
tras  850 
tras  860 
tras  870 
tras  880 
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IF(IDV.LT.5)IDV=5 

ENDIF 

IWRITE=IV1+IWIDM1 

IWIDTH=IWIDM1+1 

C 

C  PERFORM  TRIANGULAR  SLIT  INITIALIZATION.  TRANSMITTANCES  AT  A 

C  GIVEN  FREQUENCY  CONTRIBUTE  TO  2*IWIDTH-1  TRIANGULAR  SLITS. 

C  THESE  CONTRIBUTIONS  ARE  STORED  IN  ARRAY  SLIT.  WGT  IS  THE 

C  NORMALIZED  WEIGHT  USED  TO  DEFINE  THE  TRIANGLE. 

NWGT=2*IWIDTH 

WNORM=l.  /  (IWIDTH*IWIDTH) 

DO  20  I=1,IWIDTH 
WGT(I)=I*WNORM 
20  WGT(NWGT-I)=wgt(i) 

NWGT=NWGT-1 

NWGTM1=NWGT-1 

c 

c  Initialize  ICALL  (=  0  for  initial  call  to  subroutine  source) 

ICALL=0 
c 

c  Initialize  transm  (.true,  for  transmittance  only  calculations) 

transin= .  true . 

IF(IEMSCT.EQ. 1  .OR.  lEMSCT. EQ. 2) transm=. false, 
c 

c  Print  headers 

I F ( lEMSCT . EQ . 0 ) THEN 

WRITE (IPU, ’ (46H  \FREQ  TOTAL  H20  C02+  OZONE  TRACE 

1  49H  N2  CON  H20  CON  MOL  SCAT  AER-HYD  HN03  AER-HYD)  *  ) 

WRITE (IP7 ,* (46H  \FREQ  TOTAL  H20  C02+  OZONE  TRACE 

1  49H  N2  CON  H20  CON  MOL  SCAT  AER-HYD  HN03  AER-HYD)  *) 

WRITE  (IPRl,  *  (45H  \FREQ  H20  03  C02  CO  CH4 

1  47H  N20  02  NH3  NO  N02  S02 ,  / 

2  55H  \1/CM  TRANS  TRANS  TRANS  TRANS  TRANS  TRANS, 

3  39H  TRANS  TRANS  TRANS  TRANS  TRANS)*) 

WRITE(IP8,  *  (45H  \FREQ  H20  03  C02  CO  CH4 

1  47H  N20  02  NH3  NO  N02  S02 ,  / 

2  55H  \1/CM  TRANS  TRANS  TRANS  TRANS  TRANS  TRANS, 

3  39H  TRANS  TRANS  TRANS  TRANS  TRANS)*) 

ELSE  IF (lEMSCT  .EQ.  3)  THEN 

WRITE  (IPU, *(32H  \FREQ  TRANS  SOL  TR  SOLAR)*) 

WRITE  ( IP7 , *  ( 2H  \ f  T2 5 , 

+  '  40HIR:W^DIANCE  (W  M-2)  PASSED  THROUGH  FILTER, 

+  12)*)  IFILTER 

WRITE  (IP7,*(32H  \FREQ  TRANS  SOL  TR  SOLAR)*) 

ELSE  IF  (lEMSCT  .EQ.  1)  THEN 

WRITE(IPU, * (45H  \FREQ  TRANS  ATMOS.  RAD.,  T88, 

+  18H-  LOG  TOTAL  TRANS . ) * ) 

WRITE(IP7, * (2H  \,  T25, 

+  43HRADIANCE  (W  M-2  SR-1)  PASSED  THROUGH  FILTER, 

+  12)*)  IFILTER 

WRITE(IP7, * (30H  \FREQ  TRANS  ATMOS.  RAD.,  T88, 

+  18H-  LOG  TOTAL  TRANS . ) * ) 

ELSE  IF  (lEMSCT  .EQ.  2)  THEN 

WRITE (IPU, * (42H  \FREQ  TRANS  ATMOS  PATH  SINGLE, 

1  28H  GROUND  DIRECT  TOTAL  RAD)*) 

WRITE(IP7, * (2H  \,  T25, 

+  43HRADIANCE  (W  M-2  SR-1)  PASSED  THROUGH  FILTER, 
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tras  890 
tras  900 
tras  910 
tras  920 
tras  930 
tras  940 
tras  950 
tras  960 
tras  970 
tras  980 
tras  990 
traslOOO 
traslOlO 
trasl020 
trasl030 
trasl040 
traslOSO 
trasioeo 
trasl070 
trasioeo 
trasl090 
trasllOO 
traslllO 
trasll20 
trasll30 
trasll40 
trasll50 
traslieo 


trasll70 

trasllSO 

trasll90 


trasl2l0 

trasl220 


trasl230 


trasl240 

trasl250 


on  n  o  o  o  o 


+  12)*)  IFILTER 

WRITE ( IP7, * (42H  \FREQ  TRANS  ATMOS  PATH  SINGLE, 

1  28H  GROUND  DIRECT  TOTAL  RAD,  T88, 

+  18H-  LOG  TOTAL  TRANS.) ') 

END  IF 

If  (PathCard)  then 

WRITE (IPath, * (1H\,  T25, 

+  52HRADIANCE  (W  M-2  SR-1  (CM-l)--l)  PASSED  THROUGH  FILTER, 

+  I2,/,1H\,/,40H\  V  T  f  Npath*f  INTEGRAL, 

+  /,1H\)*)  If liter 

Else  if  (LastSky)  then 

WRITE { Isky, » (1H\,  T25, 

+  52HRADIANCE  (W  M-2  SR-1  (CM-l)-l) 

+  I2,/,1H\,/,40H\  V  T  f 

+  22H  Nsea*T*f  INTEGRAL, /, 1H\ )* )  Ifilter 

Else  if  (SunCard)  then 

WRITE ( Isun, » (1H\,  T25, 

+  52HRADIANCE  (W  M-2  SR-1  (CM-l)-l)  PASSED  THROUGH  FILTER, 


PASSED  THROUGH  FILTER, 
Nsky*T*f  INTEGRAL, 


I2,/,1H\,/,40H\  V 
/,1H\)*)  Ifilter 

End  If 

IF (NOPRNT . EQ . -1 ) THEN 
IF(IMULT.EQ.1)THEN 
WRITE (IPRl, * (37H 
50H  DFLX 

WRITE (IPS,  '(37H 
5 OH  DFLX 

ELSE 


Nsun*T*f  INTEGRAL, 


\V  ALTl 
DFLXS 
\V  ALTl 
DFLXS 


UFLX 

DIRS 

UFLX 

DIRS 


IF ( lEMSCT . GT . 0 ) WRITE (IPRl , 

*  (23H 

\V  ALTl 

1 

30H  B(V,T)  DTAU 

TAU) 

') 

WRITE (IPS, 

*  (23H 

\V  ALTl 

1 

30H  B(V,T)  DTAU 

TAU) 

■) 

UFLXS , 
TRANS) 
UFLXS, 
TRANS) 

ALT2, 

ALT2, 


ENDIF 

ENDIF 

Initialize  layer  loop  variables 
loopO=.true. 

call  loop(loopO,  iv,  ivx,  ikinx,inxfreq,  sumitis, transm,  iph, 

1  sumssr,  ivtest,unif, trace,  sumv,  isourc,  iday,  anglem,  frac) 
loopO=. false . 


END  INITIALIZATION,  BEGIN  OF  FREQUENCY  LOOP 

*'IVX**  IS  THE  FREQUENCY  AT  WHICH  TRANSMITTANCE  WILL  BE  CALCULATED 
DURING  THE  FIRST  PASS,  **IVX»  AND  *’IV"  MUST  BE  EQUAL. 

30  IVX=IVX+IDVX 

IF ( IV. LT. IVX) THEN 
IV=IV+IDV5 
IVTEST=.TRUE. 

ELSE 

IVTEST=- FALSE. 

ENDIF 

SET  INTERPOLATION  FRACTION. 

FRAC-FLOAT ( IV-IVX) /IDV5 
IF ( ICOUNT . EQ . ipr int ) THEN 


trasl270 
trasl280 
trasl290 
) trasl300 


trasl310 

trasl320 

trasl330 


trasl340 

trasl350 

trasl360 

trasl370 

trasl380 

trasl390 

trasl400 

trasl410 

trasl420 

trasl430 

trasl440 

trasl450 

trasl460 

trasl470 

trasl480 

trasl490 

traslSOO 

traslSlO 

trasl520 

trasl530 

trasl540 

traslSSO 

trasl560 

trasl570 

traslSSO 
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C 


1 

2 

3 

4 

5 

6 

7 

8 
9 


c 

c 

c 

c 


c 

C 


c 

C 


Reinitialize  counter  and  print  header 
ICOUNT=0 

IF { lEMSCT . EQ . 0 ) THEN 

WRITE ( IPR, * ( IHl , / 33H  FREQ  WAVELENGTH 


TOTAL 


trasl590 
trasl600 
trasl610 
H2O,trasl620 


47H  C02+ 

47H  MOL  SCAT 
//43H  1/CM 
44H  TRANS 
4 OH  TRANS 


OZONE 
AER-HYD 
MICRONS 
TRANS 
TRANS 


TRACE  N2  CONT  H20  CONT,  trasl630 
HN03  AER-HYD  INTEGRATED,  trasl640 


TRANS 

TRANS 

ABS 


TRANS  TRANS, 
TRANS  TRANS , 
ABSORPTION,/) * ) 


trasl650 
trasl660 
trasl670 
trasl680 
trasl690 
trasl700 
trasl710 
trasl720 
trasl730 
trasl740 
trasl750 
trasl760 
trasl770 
trasl780 
trasl790 
traslSOO 
traslSlO 
trasl820 
trasl830 
trasl840 
trasl850 
trasl860 
trasl870 
trasl880 
trasi890 
trasl900 
trasl910 
trasl920 
trasl930 

25HGROUND  REFLECTED  RADIANCE, TlOO , 14HTOTAL  RADIANCE, tras 19 40 
T118 , 8HINTEGRAL, T127 , 5HT0TAL, /T45 , SHTOf f  L,  T59 ,  trasl950 

6HS  SCAT,T75,5HT0TAL,T89,6HDIRECT,/1X,6H(CM-1)  ,T9,  trasl960 
7H(MICRN)  ,T19,6H(CM-1) ,T29 , 7H (MICRN)  , T39 , 6H (CM-1) ,  trasl970 
T49,7H(MICRN) , T59 , 6H (CM-1) ,T69 , 6H (CM-1) ,T79,  trasl980 

7H  (MICRN)  ,T89,6H(CM-1)  ,T99 , 6H  (CM-1)  ,  T109 , 7H  (MICRN)  ,  trasl990 
T119,6H(CM-1) ,T127,5HTRANS,/) »)  tras2000 

tras2010 
tras2020 
tras2030 
tras2040 
tras2050 
tras2060 

For  transmission  calculations,  skip  over  layer  loop  in  tratras2070 
IKMAX=1  tras2080 

ELSEIF(IMULT-EQ.l  .and.  .not.  lsame)THEN  tras2090 

tras2100 

FOR  MULTIPLE  SCATTERING  SET  IKMAX  TO  IMSMX  tras2110 

IKMAX=IMSMX  tras2 12  0 

ELSE  tras2130 

tras2140 

IF  NOT  MULTIPLE  SCATTERING,  RESET  IKMAX  TO  ORIGINAL  VALUE  tras2150 


ELSEIF ( lEMSCT . EQ . 1 ) THEN 

WRITE(IPR,  *  (1H1,20X,28HRADIANCE(WATTS/CM2-STER-XXX)  , 
/8H0  FREQ,T10,6HWAVLEN,T19,14HATMOS  RADIANCE, 

T3  9 , 9H  INTEGRAL , T4  9 , 5HTOTAL , / 2X , 6H ( CM-1 )  , 

T10,7H (MICRN) , T19 , 6H (CM-1) , T29 , 7H (MICRN) , 
T39,6H(CM-1) ,T49,5HTRANS,/)  *) 

ELSEIF (lEMSCT . EQ -  3 ) THEN 

WRITE(IPR,  *  (1H1,22X,27HIRRADIANCE  (WATTS/ CM2 -XXXX)  , 

/  7H0  FREQ ,  T1 1 ,  6HWAVLEN ,  T2  3  ,  IIHTRANSMITTED , 

T4  5 , 5HSOLAR ,  T6 1 ,  lOHINTEGRATED ,  T8  0 , 5HTOTAL , 
/2X,6H(CM-1) ,T10,7H(MICRN) , T20 , 6H (CM-1) , 
T30,7H(MICRN) , T40 , 6H (CM-1) , T50 , 7H (MICRN) , 
T60,6HTRANS. , T70 , 5HSOLAR, T80 , 5HTRANS) ») 

ELSEIF ( IMULT . EQ . 0 ) THEN 

WRITE(IPR,  *  (1H1,45X,28HRADIANCE(WATTS/CM2-STER-XXX)  , 
/ 7H0  FREQ ,  T1 1 , 6HWAVLEN ,  T2 1 , 14HATMOS  RADIANCE , 

T4 1 , 14HPATH  SCATTERED ,  T61 , 16HGROUND  REFLECTED , 

T8  5 , 5HTOTAL , T9  9 , 8HINTEGRAL , Til 0 , 5HTOTAL , 
/2X,6H(CM-1) ,T10,7H(MICRN) , T20 , 6H (CM-1) , 
T30,7H(MICRN)  , T40 , 6H (CM-1)  , T50 ,.7H  (MICRN)  , 
T60,6H(CM-1) ,T70,7H(MICRN) ,T80, 6H (CM-1) , 

T90,7H (MICRN) , TlOO , 6H (CM-1) , TllO, 5HTRANS , / ) * ) 

ELSE 

WRITE(IPR,  '  (1H1,45X,28HRADIANCE(WATTS/CM2-STER-XXX)  , 
//6H0  FREQ,T10,6HWAVLEN,T20,14HATMOS  RADIANCE, T4  0, 
4HPATH,19H  SCATTERED  RADIANCE, T69 , 


ENDIF 
ENDIF 

Determine  layer  loop  maximum 
IF(transm)THEN 
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c 

c 


40 


c 

C 

C 

C 


50 

60 

C 

C 


C 

C 


70 


IKMAX==IKMX 

ENDIF 

SUMV=0. 


tras2160 

tras2170 

tras2180 

tras2190 


Initialize  transmission  array  tras2200 
TX<1)=1.  tras2210 
TX{2)=1,  tras2220 
TX(3)=1,  tras2230 
DO  40  K=4,KMAX  tras2240 
TX(K)=0.  tras2250 
call  loop ( loop 0,iv,ivx,i)cmx,mxfreq,suinms,transm,  tras2260 


iph, sumssr , ivtest , unit, trace, sumv, isourc, iday, anglem, frac)  tras2270 

tras2280 
tras2290 
tras2300 
tras2310 
tras2320 
tras2330 
tras2340 
tras2350 
tras2360 
tras2370 
tras2380 
tras2390 
tras2400 
tras2410 
tras2420 
tras2430 
tras2440 
tras2450 
tras2460 
tras2470 
tras2480 
tras2490 
tras2500 
tran2510 

RENORMALIZE  IF  TRIANGULAR  SLIT  EXTENDS  TO  NEGATIVE  FREQUENCIEStras2520 


IF(IVX.LT.NWGTM1)THEN  tras2530 

store=l . - . 5* (NWGTMl-IVX) * (NWGTMl-IVX+1) *WNORM  tras2540 

DO  70  K=2,56  tras2550 

TX(K)=TX(K) /store  tras2560 

ENDIF  tras2570 


UNIF=TX(2) 

TRACE=TX(3) 

SDMV=TX(8) 

SDMSSR=TX(12) 

SUMMS=TX(13) 

TEB1=TX(14) 

V=FLOAT ( IVX-IWIDMl ) 
ALAM=1.0E+04/ (V+. 000001) 

Istorc  =  Istore  +  1 
Width  =  IDV*FACTOR 
f  =  Filter (V, If ilter) 

SUMA=SUMA+ (1.0 -TX ( 9 ) ) *  f *Width 
ALTX9=BIGNUM 

IF(TX(9)  .GT.O. )ALTX9=-LOG(TX(9)  ) 
GOTO (80, 90, 90, 100) ,IEMSCT+1 


tras2580 

tras2590 

tras2600 

tras2610 

tras2620 

tras2630 

tras2640 

tras2650 


tras2670 

tras2680 

tras2690 


THE  PARAMETERS  ”UNIF",  "TRACE”,  "SUMV",  "SUMSSR",  "SUMMS" 
AND  "TEBl"  ARE  TEMPORARILY  STORED  IN  "TX"  SO  THAT  THEIR 
CONVOLUTION  OVER  THE  TRIANGULAR  SLIT  CAN  BE  CALCULATED. 
TX“(2)=UNIF 
TX(3)=TRACE 
TX ( 8 ) =SUMV 
TX( 12) -SUMSSR 
TX(13)=SUMMS 
TX(14)=TEB1 
DO  60  K=2,56 
IP1=NWGT 

DO  50  I=NWGTM1, 1,-1 

SLIT (K , IPl ) =SLIT (K , I ) +WGT ( IPl ) *tx ( k ) 

IP1=I 

SLIT (K , 1 ) =WGT ( 1 ) * tx ( k ) 

TX (K) =SLIT (K, NWGT) 

CHECK  IF  VALUES  ARE  TO  BE  PRINTED 
IF ( IVX . LT . IWRITE) GOT03  0 
IWRITE=IWRITE+IDV 
IF  (IWRITE .  GT.  IVXMAX)  FACTOR= .  5 
I COUNT=I COUNT+ 1 
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c 

c 

80 


1 

1 

1 

1 

1 


c 

c 


1 

1 

1 

1 

1 


c 

c 

c 

c 


tras2700 


TRANSMITTANCE  ONLY  tras2710 

TX(10)=1.-TX(10)  tras2720 

TX(7)=TX(7)*TX(16)  tras2730 

WRITE (IPR,  * (F8.0,F8.3,11F9.4,F12.3)  * ) V, ALAM, TX (9) ,TX(17) ,  tras2740 

UNIF,TX(31) , TRACE, TX(4) ,TX(5) ,TX(6) ,TX(7) ,TX(11) ,TX(10) , SUMAtras2750 
WRITE (IPRl, * (F7,0,11F8.4,1PE10.3) »)V,TX(17) ,TX(31) ,TX(36) ,  tras2760 

TX(44) ,TX(46) ,TX(47) ,TX(50) ,TX(52) ,TX{54) ,TX(55) ,TX(56)  tras2770 

WRITE(IP8,  • (F7.0,11F8.4,1PE10.3) »)V,TX{17) ,TX(31) , TX(36) , 

TX(44) ,TX(46) ,TX(47) ,TX(50) ,TX(52) ,TX(54) ,TX{55) ,TX(56) 

WRITE (IPU, * (F7,0,11F8.4,1PE10.3) * ) V,TX(9) ,TX (17) ,UNIF, 

TX(31) , TRACE, TX(4) ,TX(5) ,TX(6) ,TX(7) ,TX(11) ,TX(10) ,ALTX9 
WRITE(IP7, * (F7.0,11F8.4,1PE10,3) *)V,TX(9) ,TX(17) ,UNIF, 

TX(31) , TRACE, TX{4) ,TX(5) ,TX(6) ,TX(7) ,TX(11) ,TX(10) ,ALTX9 
GOTOllO 


tras2780 

tras2790 


tras2800 

tras2810 


ATMOSPHERIC  RADIANCE  INCLUDING  EMISSION  OF  BOUNDARY 
ATTENUATED  BY  TOTAL  TRANSMISSION 


tras2820 

tras2830 

tras2840 


CALCULATE  THERMAL  RADIANCE  CONTRIBOTION  OF  BOUNDARY  AND  tras2850 

ADD  THE  SCATTERED  CONTRIBUTION  TO  THE  THERMAL  RADIANCE  tras2860 

IF  THE  PATH  INTERSECTS  THE  SURFACE  tras2870 

IF  ( (TBOUND.LE.O. )  .OR.  (PathCard) )  THEN 

BBG  =  0.  tras2890 

ELSE  tras2900 

BBG=BBFN(TBOUND,V) *TX(9) *EMISS  tras2910 

IF(IMULT.EQ.l  .AND.  ground)  THEN 

BBG=BBG+SALB*FDNTRT*TX(9) /PI  tras2920 

END  IF 

ENDIF  tras2930 

tras2940 


ADD  THERMAL  BOUNDARY  AND  MULTIPLE  SCATTERED  RADIANCE 
SUMV= (SUMV+BBG) *f . 

SUMW=SUMV 
IF  (V.GT.O.)  THEN 

SUMV=(1.0E+08/V**2)*SUMV  1  W  m-2  sr-1  (cm-l)-l 

END  IF 

IF ( lEMSCT . EQ . 1 ) THEN 

RADSUM=RADSUM  +  SUMV*Width 

WRITE (IPR, » (F8,0,F8.3,1P3E10.2,0PF9.4) *) 

V ,  ALAM ,  SUMV ,  S  UMW ,  RA  D  S  UM ,  T  X  ( 9 ) 

WRITE (IPU, * (F7.0,F8.4,1PE15.8,T96,E10.3) ’) 

V,TX(9) ,SUMV,ALTX9 

WRITE (IP7, * (F7.0,F8.4,1PE15.8,T96,E10.3) ') 

V,TX(9) ,SUMV,ALTX9 

WRITE (IPrl, ’ (F7.0,11F8.4,1PE10.3) ‘ ) V,TX (9) , TX ( 17 ) ,UNIF, 

TX(31) , TRACE, TX(4) ,TX(5) ,TX(6) ,TX(7) ,TX(11) ,TX(10) ,ALTX9 
WRITE (IP8,  * (F7.0,11F8,4,1PE10,3) * ) V,TX(9) ,TX(17) ,UNIF, 

TX(31) , TRACE, TX(4) ,TX(5) ,TX(6) ,TX(7) ,TX(11) ,TX(10) ,ALTX9 
SUMT=SUMV 
SUMTT=SUMW 

ELSE 


tras2950 

tras2970 


tras2990 

tras3010 

tras3020 

tras3030 

tras3040 


tras3050 

tras3060 


tras3070 

tras3080 

tras3090 

tras3100 


SOLAR  SCATTERED  RADIANCE  tras3110 

CALL  SOURCE (V, ISOTOC, IDAY , ANGLEM, SS)  tras3120 

tras3130 


MULTIPLY  SUMSSR  BY  THE  EXTRATERRESTRIAL  SOURCE  STRENGTH  SStras3140 
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SUMSSS=SUMSSR*SS 


CALCULATE  TOTAL  SINGLE  SCATTERED  +  MULTIPLE  SCATTERED 


tras3150 

tras3160 

tras3170 


SOLAR  RADIANCE  FOR  EACH  FREQUENCY  [W/CM2-STER-MICROMETER]  tras3180 
SUMSSR=SUMSSS+SUMMS  tras3190 

store=0.  tras3200 

if  (v,gt.O,  )store=l.e8/v**2  1  [W  iti~2  sr-1  (cra-l)-l] 

SUMS=store*SUMSSR  tras3220 

SUMSSS=store*SUMSSS  tras3230 


RFLSOL  IS  GROUND-REFLECTED  DIRECT  SOURCE  RADIANCE  AND 
RFLS0L=0 . 

RFLS=0 . 

RFLSS=0. 

RFLSSS=0. 

IF(ground  .AND.  TEBl .GT. 0) THEN 


tras3l90 

tras3200 

tras3220 

tras3230 

tras3240 

tras3250 

tras3260 

tras3270 

tras3280 

tras3290 

tras3300 


tras3400 

tras3410 

tras3420 

tras3430 

tras3440 

tras3450 

tras3460 

tras3470 

tras3480 

tras3490 


IF  ( ANGSUN .  GE .  0 .  )  RFLSSS=SS*TEB1*SALB*C0S  ( ANGSUN*CA)  /PI  tras3 310 
RFLSOL=RFLSSS  tras3320 

IF(IMULT.EQ.  l)RFLSOL=RFLSOL+SALB*FDNSRT*TX{9)  /PI  tras3330 

RFLS=STORE  *RFLSOL  tr as  3340 

RFLSS=STORE*RFLSSS  tras3350 

ENDIF  tras3360 

SUMT=SUMV+ ( SUMS+RFLS ) *  f 
SUMTT=SUMW+(SUMSSR+RFLSOL)  *f 
RADSUM=RADSUM  +  SUMT*Width 

IF  (IMULT.NE.l)  THEN  tras3400 

WRITE (IPR, ’ (F8.0,F8.3,1P9E10.2,0PF9.4) *)  tras3410 

V ,  ALAM ,  SUMV ,  SUMW ,  SUMS ,  SUMSSR ,  RFLS ,  RFLSOL ,  tr  as  3  420 

SUMT ,  SUMTT ,  RADSUM ,  TX  ( 9 )  t r as  3 4  30 

else  tras3440 

WRITE (IPR, * (F7.0,F8.3,1P11E10.2,0PF7,4) *)  tras3450 

V ,  ALAM ,  SUMV ,  SUMW ,  SUMS ,  SUMS  SR ,  SUMS  S  S ,  RFLS ,  tr  a  S  3  4  6  0 

RFLSOL , RFLSS , SUMT , SUMTT , RADSUM , TX ( 9 )  tr as  3  470 

END  IF  tras3480 

WRITE (IPU,  » (F7.0,F8.4,1P6E9.2,0P2F8.4,T96,1PE10.3)  * )V,  tras3490 

TX  ( 9 )  ,  SUMV ,  SUMS ,  SUMSSS ,  RFLS , RFLSS ,  SUMT ,  TEBl ,  TEB2SV ,  ALTX9tras3 500 
WRITE (IP7^  * (F7.0,F8.4,1P6E9.2,0P2F8.4,T96,1PE10.3)  ’)V, 

TX  ( 9 )  ,  SUMV ,  SUMS ,  SUMSSS ,  RFLS ,  RFLSS ,  SUMT ,  TEBl ,  TEB2 SV ,  ALTX9 
WRITE(IPrl, * (F7.0,11F8.4,1PE10.3) ^V^TXCO) ,TX(17) ,UNIF,  tras3510 

TX(31)  , TRACE, TX(4)  ,TX(5)  ,TX(6)  ,TX(7)  ,TX(11)  ,TX(10)  ,ALTX9  tras3520 

WRITE (IP8,  * {F7.0,11F8.4,1PE10.3) ’ ) V,TX(9) ,TX(17) ,UNIF, 

TX(31) , TRACE, TX(4) ,TX(5) ,TX(6) ,TX(7) ,TX(11) ,TX(10) ,ALTX9 
END  IF  tras3530 

IF  (PathCard)  THEN 

Tau(Istore)  =  TX(9) 

PathSum  =  PathSum  +  SUMT*Width 
Write(Ipath,  > (I5,2F7.3,2 (1PE11.3) ) *) 

V,  TX(9),  f,  SUMT,  PathSum 

END  IF 

IF  (SkyCard)  skyN{Msea, Istore)  =  SUMT 

GO  TO  110  tras3540 

tras3550 

DIRECTLY  TRANSMITTED  SOLAR  IRRADIANCE  [WATTS/ (CM2  MICROMETER)  ]tras3 5 60 
CALL  SOURCE  (V,  ISOURC,  IDAY ,  ANGLEM,  SOLIL)  tras3570 

SOLIV=0.  tras3580 

IF(V.GT.O.)SOLIV=SOLIL*l.E+8/V**2  I  [W  m-2  sr-1  (cm-l)-l] 
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TS0LIV=S0LIV*TX(9)  *f 
IF  (SunCard)  then 

Ho(lstore)  =  TSOLIV 
END  IF 

TSOLIL=SOLIL*TX { 9 )  *f 

STSOL=STSOL+TSOLIV*Width 

SSOL=SSOL+SOLIV*Width 

WRITE(IPR, * (F8.0,F8.3,1P6E10.2,0PF9,4) ») 

1  V ,  ALAM ,  TSOLIV ,  TSOLIL ,  SOLI V ,  SOLIL ,  STSOL ,  SSOL ,  TX  ( 9  ) 

WRITE (IPU,  '  (F7,0,F8,4,1P2E9.2,T96,E10.3)  *) 

1  V ,  TX  ( 9 )  ,  TSOLIV ,  SOLIV ,  ALTX9 

WRITE(IP7; * (F7.0,F8,4,1P2E9.2,T96,E10.3) ») 

1  V , TX ( 9 ) , TSOLIV , SOLIV , ALTX9 

WRITE (IPrl,  '  (F7.0,11F8.4,1PE10.3)  *)V,TX(9) ,TX(17) ,UNIF, 

1  TX(31)  , TRACE, TX(4) ,TX(5) ,TX(6)  ,TX(7) ,TX(11) ,TX(10) ,ALTX9 

WRITE (IP8,  * (F7,0,11F8.4,1PE10.3) * ) V,TX(9) ,TX (17) ,UNIF, 

1  TX(31)  , TRACE, TX(4)  ,TX(5)  ,TX(6)  ,TX(7)  ,TX(11)  ,TX(10)  ,ALTX9 

SUMT=TSOLIV 
RADSUM=STSOL 

110  IF ( lEMSCT - NE . 0 ) THEN 

IF  (SUMT .  GE .  RADMAX)  THEN 
VRMAX=V 
RADMAX=SUMT 
ENDIF 

IF  (SUMT .  LE  -  RADMIN)  THEN 
VRMIN=V 
RADMIN==SUMT 
ENDIF 
ENDIF 
FACTOR=l. 

IF  ( IWRITE .  LE .  IVXMAX)  G0T03  0 

END  OF  FREQUENCY  LOOP 

IVX=INT (V+ . 5 ) 

IF  (IFILTER  -EQ.  0)  THEN 
Sumf  =  IVX  -  IVl 

ELSE  IF  ((1  .LE.  IFILTER)  .AND.  (IFILTER  .LE.  6))  THEN 
Sumf  =  FLIST(5,  IFILTER) 

END  IF 


tras3640 

tras3650 

tras3660 

tras3670 


tras3680 

tras3690 


tras3700 

tras3710 

tras3720 

tras3730 

tras3740 

tras3750 

tras3760 

tras3770 

tras3780 

tras3790 

tras3800 

tras3810 

tras3820 

tras3830 

tras3840 

tras3850 

tras3860 

tras3870 


IF  ((.NOT.  Sea)  .AND.  ({lEMSCT  .EQ.  1)  .OR.  (lEMSCT  .EQ.  2)))  THEN 
IF  (IFILTER  ,EQ.  0)  THEN 
VI  =  FLOAT  (IVl) 

V2  =  FLOAT  (IVX) 
dV  =  FLOAT  (IDVX) 

CALL  RtoT  (VI,  V2,  dV,  RADSUM/1.E4,  BBTEMP) 

ELSE  IF  ((IFILTER  .GE.  1)  .AND.  (IFILTER  .LE.  6))  THEN 
CALL  Planck  (RADSUM/1, E4 ,  IFILTER,  BBTEMP) 

END  IF 
END  IF 


IF  (Sea)  THEN 

IF  (PathCard)  THEN 

Npath  =  PathSum 
PathTrans  =  1. -SUMA/Sumf 
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PathRange  =  RANGE 
PathAngle  =  ANGLE 
ELSE  IF  (LastSky)  THEN 
DO  I  =  1,  Istore 
DO  K  =  1,  3 

Rsky(K)  =  SkyN(K,I) 

END  DO 

CALL  Fit (Tsky,Rsky,3,a,b) 

V  =  IVl  +  (I  -  l.)*IDV 
f  =  FILTER (V, If liter) 

IF  ((I  .EQ,  1)  .OR.  (I  .EQ.  Istore))  THEN 
Width  =  IDV/2. 

ELSE 

Width  =  IDV 

END  IF 

Nbb  =  BBFN(TBOUNDold,V)*f*lE8/V**2 
CALL  Sky ( Tr , Pr , WSS , a , b , V , Bover A , ev , Nvsky ) 

Nsky  =  Nsky  +  Nvsky*Width*Tau(I) 

Nsea  =  Nsea  +  ev*Nbb*Width*Tau(I) 

Write(Isky,  ’(15,  2F7,3,  4 (IPEll. 3) ) ’ ) 

+  V,  Tau(I) ,  f,  Nvsky*Tau(I) ,  Nsky, 

+  ev*Nbb*Tau (I) ,  Nsea 

END  DO 

ELSE  IF  (SunCard)  THEN 
DO  I  =  1,  Istore 

V  =  IVl  +  (I  -  l.)*IDV 
f  =  FILTER (V,  If liter) 

IF  ((I  -EQ.  1)  .OR.  (I  .EQ.  Istore))  THEN 
Width  =  IDV/2, 

ELSE 

Width  =  IDV 

END  IF 

CALL  Sun(Tr,Pr,WSS,To,Po,V,suin4) 

No  =  Ho (I) *sum4/ (4 . *pi*epsilon**2*Bov€rA) 

Nsun  =  Nsun  +  No*Width*Tau (I) 

Write (Isun,  ’ (I5,2F7.3,2(1PE11.3) ) *) 

+  V,  Tau(I),  f,  No*Tau(I) ,  Nsun 

END  DO 

END  IF 

Ntotal  =  Npath  +  Nsea  +  Nsky  +  Nsun 
IF  (IFILTER  .EQ.  0)  THEN 
VI  =  FLOAT  (IVl) 

V2  =  FLOAT  (IVX) 
dV  =  FLOAT  (IDVX) 

CALL  RtoT  (VI,  V2,  dV,  Ntotal/ 1.E4,  TotalT) 

ELSE  IF  ((IFILTER  .GE,  1)  .AND.  (IFILTER  .LE.  6))  THEN 
CALL  Planck  (Ntotal/l.E4 ,  IFILTER,  TotalT) 

END  IF 
END  IF 
C 

IF  (.NOT.  Sea)  THEN 

WRITE (IPR, • (26HINTEGRATED  ABSORPTION  FROM, 15, 3H  TO, 15, 

+  7H  CM-1  =,F10.2,5H  CM-1 , /23HAVERAGE  TRANSMITTANCE  =, 

+  F6,4,/)*)  IVl,IVX,SUMA,l.-SUMA/SuTnf 

WRITE (IP6,  * ( 

+  /24HINTEGRATED  ABSORPTION  =,  F10.2, 

+  lOH  CM-1  FROM,  15,  3H  TO,  15,  5H  CM-1, 


tras3880 

tras3890 
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+ 

+ 

+ 

+ 

+ 
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+ 

+ 

+ 

+ 

+ 
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+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 


+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 


/ 2 4HAVERAGE  TRANSMITTANCE  =,  F12.4)*) 

SUMA ,  I VI ,  I VX ,  1 .  -SUMA/ Suinf 
IF  (lEMSCT  .EQ.  0)  THEN 

WRITE (IP4,  * (T13,F7.2,F8.3,F10,3,F8.3)  *) 

(90. -ANGLE) *pi/180.*l.E3, 

ANGLE ,  RANGE ,  1 . -SUMA/ Sumf 
ELSE  IF  ((lEMSCT  .EQ.  1)  .OR.  (lEMSCT  .EQ.  2))  THEN 
WRITE (IPR, * (22  HINTEGRATED  RADIANCE  =,1PE11.3, 
lOH  WATTS  M-2,7H  STER-1, 

/22H  MINIMUM  RADIANCE  =,E11.3,10H  WATTS  M-2 , 
19H  STER-1  (CM-l)-l  AT, OPFll . 1, 5H  CM-1, 

/8H  MAXIMUM, 14H  RADIANCE  =,1PE11.3, 

29H  WATTS  M-2  STER-1  (CM-l)-l  AT , OPFll . 1 , 5H  CM-1 
23H  BOUNDARY  TEMPERATURE  =,  F11.2,2H  K, 

/22H  BOUNDARY  EMISSIVITY  =,F12.3)  *) 

RADSUM ,  RADMIN ,  VRMIN ,  RADMAX ,  VRMAX ,  TBOUND ,  EMI  SS 
WRITE(IP6,*( 

/24HMAXIMUM  RADIANCE 
23H  W  M-2  SR-1  (CM-l)-l  AT, 

/24HMINIMUM  RADIANCE 
23H  W  M-2  SR-1  (CM-l)-l  AT, 

/24HBOUNDARY  TEMPERATURE 
/24HBOUNDARY  EMISSIVITY 
/24HFILTERED  RADIANCE 
IIH  W  M-2  SR-1, 

/24HBLACKBODY  TEMPERATURE 
RADMAX ,  VRMAX ,  RADMIN ,  VRMIN ,  TBOUND ,  EMI  S S  , 

RADSUM,  BBTEMP  -  TO 

WRITE ( IP4 , * (T13 , F7 . 2 , F8 . 3 , FIO . 3 , F8 . 3 , 5 ( IPEIO . 3 ) , 

0PF8.1) * 

(90. -ANGLE) *pi/180.*l.E3,  ANGLE,  RANGE, 
l.-SUMA/Sumf ,  RADSUM,  Nsea,  Nsky,  Nsun, 

RADSUM,  BBTEMP-TO 


1PE11.3, 

0PF8.1,  5H  CM-1, 
1PE11.3, 

0PF8.1,  5H  CM-1, 
F11.2,  2H  K, 
F12.3, 

1PE11.3, 


=,  OPFll. 1,  2H  C) *) 


ELSE 


END 


END  IF 
IF 


WRITE (IPR, * (24H  INTEGRATED  IRRADIANCE  =,1PE11.3, 

lOH  WATTS  M-2,724H  MINIMUM  IRRADIANCE  =,E11.3 
13H  WATTS  M-2  AT, OPFll. 1, 5H  CM-1,/10H  MAXIMUM  , 
14H  IRRADIANCE  =,1PE11.3, 

23H  WATTS  M-2  (CM-l)-l  AT, OPFll . 1 , 5H  CM-1)*) 
RADSUM ,  RADMIN ,  VRMIN ,  RADMAX ,  VRMAX 
WRITE (IP6, * ( 

24HINTEGRATED  IRRADIANCE  =,  1PE11.3,  6H  W  M-2, 
/24HMINIMUM  IRRADIANCE  =,  Ell. 3, 

18H  W  M-2  (CM-l)-l  AT,  0PF8.1,  5H  CM-1, 
/24HMAXIMUM  IRRADIANCE  =,  1PE11.3, 

18H  W  M-2  (CM-l)-l  AT,  0PF8.1,  5H  CM-1)*) 

RADSUM ,  RADMIN ,  VRMIN ,  RADMAX ,  VRMAX 


IF 

(( (LastSky) 

.AND.  (lEMSCTold 

•EQ.  1)) 

.OR.  (SunCard) 

WRITE (IP 6, 

*(/,  24HRECEIVED 

RADIANCE 

VALUES, 

//,  T10,24H 

PATH  TO  FOOTPRINT 

~ , 

FIO. 5,  IIH 

W  M-2  SR-1, 

+ 

T56,12H 

(AV.  TRANS., 

F7.4, 

IH)  , 

+ 

/,  T10,24H 

SEA  EMISSION 

r 

FIO. 5, 

IIH  W 

+ 

/,  T10,24H 

SKY  REFLECTION 

t 

FIO. 5, 

IIH  W 

tras4020 


) 

tras3920 
,  tras3930 
tras3940 

tras3960 
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+  /,  T10,24H  SUN  GLINT  =,  F10.5,  IIH  W  M-2  SR-1 

+  //,  T10,24H  TOTAL  RADIANCE  =,  F10,5,  IIH  W  M-2  SR-1 

+  /,  T10,24H  BLACK  BODY  TEMP.  =,  FlO.l,  2H  C  ) * ) 

+  Npath,  PathTrans,  Nsea,  Nsky,  Nsun,  Ntotal,  TotalT-TO 

WRITE (IP4,  * (T13,F7,2,F8.3,F10.3,F8.3,5(1PE10.3) ,0PF8.1)  *) 

(90 . “PathAngle) *pi/180 . *1. E3 ,  PathAngle,  PathRange, 
PathTrans,  Npath,  Nsea,  Nsky,  Nsun,  Ntotal,  TotalT-TO 
Npath  =  0. 

Nsea  =  0. 

Nsky  =  0 . 

Nsun  =  0. 

END  IF 

RETURN 

END 


tras4100 

tras4110 
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*******************************  MOD22 • FOR  ****************************** 
*  * 

*  New  version  of  Cox-Munk  routines  with  integration  over  sea  ^  * 

*  slopes  and  interpolation  between  three  sky  angles  for  estimation  * 

*  of  incident  sky  radiance.  * 

*  * 

*  Last  revised:  July  14,  1995.  * 

*  * 
************************************************************************ 


SUBROUTINE  Sky ( Tr , Pr , W , a , b , v , Bover A , e , N sky) 

REAL  Nsky 
CU  USES  rho 

C  Outputs : 

C  Calculates  (1)  the  normalization  factor  "BoverA"  in  the 

C  denominator  of  the  interaction  probability,  and  spectral 

C  values  for  (2)  the  effective  emmisivity  *'e”  of  the  ocean 

C  surface,  and  (3)  the  sky  radiance  "Nsky”  [W  m-2  sr-1  (cm-l)-l] 

C  reflected  from  the  ocean  surface. 

C  Inputs: 

C  The  receiver  spherical  coordinates  [rad]  are  (Tr,Pr) .  The 

C  wind  speed  is  W  [m  s-1] .  v  [cm-l]  is  the  wavenumber,  a  and  b  are 

C  coefficients  of  a  least  squares  fit  such  that  Ns,  the  spectral 

C  sky  radiance  [W  m-2  sr«l  {cm-l)-l]  incident  on  the  ocean 

C  at  zenith  angle  Ts  [rad]  ,  is  given  by 

C  Ns(Ts,v)  =  l./[a(v)  -  b(v)*Ts**2], 

Last  revision: 

January  27,  1995. 

COMMON/ Constants /pi , r2d , d2r , epsilon , delta , onem, onep , infinity 
REAL  infinity.  Ns 
if  (W  .ge.  .01)  then 

use  the  Cox-Munk  standard  deviation  for  a  real  sea 
Su  =  sqrt(3.16E-3*W) 

Sc  =  sqrt(3.E-3  +  1.92E-3*W) 

else 

use  a  delta  function  for  an  ideal  calm  sea 
Su  =  .01 
Sc  =  .01 

end  if 

pO  =  1./ (2.*pi*Su*Sc) 

Sav  =  (Su  +  Sc) /2 . 

N  =  2 
M  =7 

Smx  =  N*2.303*Sav 
dS  =  Smx/M 
Ar  =  sin(Tr) *cos(Pr) 

Br  =  sin(Tr) *sin(Pr) 

Cr  =  cos(Tr) 

suml  =  0. 
sum2  =  0. 
sum3  =  0. 


C 


C 
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do  Sx  =  -rSmx,  Smx,  dS 

Symx  =  sqrt  (abs  (Sinx**2  -  Sx**2)) 
do  Sy  =  -Symx,  Symx,  dS 

For  each  position  (Sx,Sy)  in  slope  space: 

calculate  the  occurrence  probability  density  p: 
arg  =  ((Sx/Su)**2  +  (Sy/Sc) **2) /2 . 
if  ((arg)  .ge.  log (pO/delta) )  then 

p  -  0. 

else 

p  =  pO*exp(-arg) 

end  if 

calculate  omega,  the  angle  of  incidence  and  Ts, 

the  zenith  angle  of  the  source  ray. 

dd  =  Sx**2  +  Sy**2 

fO  =  -  Ar*Sx  -  Br*Sy  +Cr 

w  =  fO/sqrt(l.  +  dd) 

if  ((onem  .le.  w)  .and.  (w  .le.  onep)  )  then 
omega  =0. 

else  if  ((-onep  .le.  w)  .and.  (w  .le.  -onem))  then 
omega  =  pi 

else 

omega  =  acos(w) 

end  if 

uu  =  (-  2.*Ar*Sx  -  2.*Br*Sy  +  Cr*(l.  -  dd))/(l.  +  dd) 
if  ((onem  .le.  uu)  .and.  (uu  .le.  onep))  then 
Ts  =  0. 

else  if  ((-onep  .le.  uu)  .and.  (uu  .le.  -onem))  then 
Ts  =  pi 

else 

Ts  =  acos(uu) 

end  if 

interpolate  for  Ns(Ts), 

Ns  ==  l./(a  -  b*(Ts**2)) 

define  integrands, 
fl  =  f0*p 

f2  =  rho (omega, v) *fl 
f3  =  Ns*f2 

and  accumulate  integrals  over  all  slopes, 
if  (omega  .le.  pi/ 2.)  then 
suml  =  fl  +  suml 
sum2  =  f2  +  sum2 

if  (Ts  .le.  pi/2.)  sum3  =  f3  +  sum3 
end  if 
end  do 
end  do 

suml  =  suml*dS**2 
sum2  =  sum2*dS**2 
sum3  =  sum3*dS**2 

BoverA  =  suml 
e  =  1.  -  sum2/suml 
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Nsky  =  suin3/suml 


return 

END 

SUBROUTINE  Sun ( Tr , Pr , W , To , Po , V , sum4 ) 

CU  USES  rho 

Outputs : 

Calculates  a  spectral  solar  reflectivity  **suin4”  for  the 
ocean  surface  apart  from  a  normalization  factor  of 
(4 . *”BoverA” )  or  (4.*”suml”). 

Inputs: 

The  receiver  spherical  coordinates  [rad]  are  (Tr,Pr) .  The 
wind  speed  is  W  [m  s-l].  The  spherical  coordinates  [rad] 
of  the  solar  center  are  (To,Po) .  v  [cm-l]  is  the  wavenumber. 

Note: 

The  larger  the  value  of  M,  the  y  coordinate  step  size,  the 
more  precise  and  slower  the  sum.  For  fixed  M  precision 
improves  with  wind  speed.  For  W  =  1  m  s-1  and  M  =  5,  the 
precision  is  better  than  1.5  %  around  the  center  of  the 
glint  pattern  until  the  receiver  zenith  angle  exceeds  89.5 
degrees . 

Bug: 

Divides  by  zero  when  the  sun  is  on  the  zenith. 

Last  Revision: 

January  27,  1995. 

COMMON/ Constants /pi , r 2d , d2r , epsilon , delta , onem , onep , inf inity 
REAL  infinity 

C  Find  the  rectangular  receiver  coordinates 
Ar  ==  sin(Tr)  *cos  (Pr) 

Br  =  sin(Tr) *sin [Pr) 

Cr  =  cos(Tr) 

C  Find  the  Cox-Munk  wind  dependent  slope  standard  deviations 
if  (W  .ge.  .01)  then 

C  use  the  Cox-Munk  standard  deviation  for  a  real  sea 

Su  =  sqrt(3.16E-3*W) 

Sc  =  sqrt(3.E~3  +  1.92E-3*W) 

else 

C  use  a  delta  function  for  an  ideal  calm  sea 

Su  =  .01 
Sc  =  .01 

end  if 

pO  =  1./ (2.*pi*Su*Sc) 

M  =  5 

N  =  2*M  +  1 
sum  =  0 . 

dY  =  2.*epsilon/N 
do  I  -  1,  N 
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Y  =  epsilon  -  (I  -  0.5) *dY 
Xmax  =  sqrt (epsilon**2  -  Y**2) 

K  =  int(2.*Xinax/dY  +  onein*0.5) 

dX  =  2 . *Xmax/K 

do  J  =  1,  K 

X  =  Xmax  -  (J  -  0.5) *dX 

For  each  position  (X,y)  (rectangular  coordinates  with 
respect  to  the  solar  center)  on  the  solar  disk, 

Find  the  spherical  source  coordinates: 

Ts  =  To  -  y 

Ps  =  Po  -  X/sin(To) 

if  (Ts  .gt.  pi/2.*onep)  then 

print  *,  "Error  from  'Sun*:  Part  of  solar  disk" 
print  "  is  below  horizon.” 

return 
endif 

Find  the  the  rectangular  source  coordinates: 

As  =  sin(Ts) *cos (Ps) 

Bs  =  sin(Ts) *sin(Ps) 

Cs  =  cos(Ts) 

Find  the  slopes  (Sx,Sy)  for  a  specular  reflection  from 
source  (Ts,Ps)  to  receiver  (Tr,Pr) : 
if  (abs(As  +  Ar)  .le.  delta)  then 
Sx  =  0. 

else  if  ((Cs  +  Cr)  -le.  delta)  then 
Sx  =  sign(inf inity ,  -(As  +  Ar) ) 

else 

Sx  =  -  (As  +  Ar)/(Cs  +  Cr) 

end  if 

if  (abs(Bs  +  Br)  .le.  delta)  then 
Sy  =  0. 

else  if  ((Cs  +  Cr)  .le.  delta)  then 
Sy  =  sign(inf inity,  -(Bs  +  Br) ) 

else 

Sy  =  -  (Bs  +  Br)/(Cs  +  Cr) 

end  if 

Find  the  Cox-Munk  occurrence  probability  density: 
arg  =  ((Sx/Su)**2  +  (Sy/Sc) **2) /2 . 
if  ((arg)  .ge.  Iog(p0/delta) )  then 

p  =  0. 

else 

p  =  p0*exp(-arg) 

end  if 

Find  the  angle  of  incidence  (omega)  for  the  specular 
reflection: 

dd  =  (1,  +  As*Ar  +  Bs*Br  +  Cs*Cr)/2. 
if  (dd  .le.  delta)  dd  =  0. 
ss  =  sqrt(dd) 

if  ((onem  .le.  ss)  .and.  (ss  .le.  onep) )  then 
omega  =  0. 

else  if  ((-onep  .le.  ss)  .and.  (ss  .le.  -onem))  then 
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(A12) 


(Ai3) 


{A14) 

(A14) 


non 


(A14) 


c 

c 


C 

c 

c 

c 


omega  =  pi 

else 

omega  =  acos(ss) 

end  if 

Find  the  facet  tilt  (Tn)  for  the  specular  reflection: 

Tn  =  atan(sqrt(Sx**2  +  Sy**2))  (A6) 

Integrate : 

sum  =  rho (omega, v) / (cos (Tn) **4) *p*dX  +  sum  (32) 

end  do 
end  do 

sum4  =  sum*dY 

return 

END 

SUBROUTINE  Fit (x,y , n, a,b) 

DIMENSION  X(*),y(*) 

Given  (x,y)  pairs  in  the  data  arrays  x(i)  and  y(i) ,  where 
1  <=  i  <=  n,  performs  a  least  sqares  fit  of  these  data  to 
the  equation 

y  =  1/ (a  -  b*x**2) 

and  returns  the  values  of  a  and  b. 

Last  revised:  March  13,  1995. 

DOUBLE  PRECISION  nn,  bb,  cc(4) 

where  cc(l:4)  =  cOl,  c21,  c20,  c40. 

nn  =  FLOAT (n) 

do  i  -  1,  n 

if  (y(i)  -eq.  0.)  then 
a  =  7.E5 
b  =  0. 
return 
end  if 
end  do 


do  i  =  1, 

4 

cc(i) 

= 

O.dO 

end  do 

do  i  =  1, 

n 

cc(l) 

= 

cc(l) 

i-/y{i) 

cc(2) 

= 

cc{2) 

+ 

x(i)**2/y(i) 

cc(3) 

= 

cc(3) 

+ 

x(i)**2 

CC  (4) 

= 

cc(4) 

+ 

x(i)**4 

end  do 

bb  =  (nn*cc(2)  -  cc(3) *cc (1) ) / (nn*cc(4)  -  cc(3)**2) 
a  =  (cc(l)  -  cc(3)*bb)/nn 
b  =  -  bb 
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END 


SUBROUTINE  Foot (ThetaO ,  PhiO,  ThetaS,  PhiS,  PsiPO,  Beta,  Psi) 

COMMON  /Constants /  Spi , Sr 2d , Sd2r , eps i Ion , de Ita , onem , onep , infinity 
COMMON  /Geometry/  To,Po,Tr,Pr 

comm  /cLd2/  I^ZE:isksN,IVULCN,ICSTL,ICLD,IVSA,VIS,W,WHH, 

+  RAINRT 

COMMON  /IFIL/  IRD,IPR,IPU,NPR,IPR1,IP6,IP7,IP8,IP4,IRDS,IP6S, 

+  ITR, Ipath, Isky , I sun 

COMMON  /Sea/  Sea,  Hit,  Msea,  TBOUNDold,  lEMSCTold 
REAL  infinity 

************************************************************************ 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

■k 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


This  routine  calculates  zenith  angles  and  azimuthal  angles  from 
the  footprint,  defined  as  the  location  where  a  hit  has  just 
occurred . 

All  the  arguments  are  inputs,  and  all  are  angles  in  degrees; 
ThetaO  and  PhiO  are  the  observer  latitude  and  longitude. 
ThetaS  and  Phis  are  the  solar  latitude  and  longitude. 

PsiPO  is  the  path  azimuth  (+  E  of  N)  seen  by  the  observer - 
Beta  is  the  angle  subtended  at  the  center  of  the  earth 
between  the  observer  and  the  footprint - 
PsiW  is  the  wind  azimuth  (+E  of  N)  seen  by  the  observer. 

The  outputs  are  To,  Po,  and  Pr,  angles  in  radians  passed 
through  the  common  block  "Geometry” ; 

To  is  the  zenith  angle  of  the  center  of  the  sun 
from  the  footprint. 

Po  is  the  azimuth  (+  W.  of  PsiW)  of  the  center  of  the  sun 
from  the  footprint. 

Pr  is  the  azimuth  (+  W.  of  PsiW)  receiver  as  seen 
from  the  footprint. 

Note:  Tr  has  been  calculated  in  DPFNMN  and  is  used 

in  this  subroutine  only  for  printing  to  "OUT". 

Last  revision:  June  15,  1995. 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

★ 

•k 

* 

* 

* 

* 

* 

* 

* 


************************************************************************ 

DOUBLE  PRECISION  DThGtaO,  DPhiO,  DThetaS,  DPhiS,  DPsiPO,  DBeta, 

+  DPsiW,  DThetaF,  DPhiF,  Pi,  Side,  Angle,  DTo,  DPo 


C 


Pi  =  4.*DATAN(1.) 

D2R  =  Pi/180. 

R2D  =  REAL (180. /Pi) 

First,  convert  to  radians  and  increase  precision: 


DThetaS  =  DBLE (ThetaS) *D2R 
DPhiS  =  DBLE ( Phis )*D2R 
DThetaO  =  DBLE (ThetaO) *D2R 
DPhiO  =  DBLE (PhiO) *D2R 
DPsiPO  =  DBLE (PsiPO) *D2R 
DBeta  =  DBLE (Beta) *D2R 
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DPsiW  =  DBLE(Psi+PsiPO) *D2R 

then  use  the  geometry  of  three  spherical  triangles  connecting 
the  north  pole,  the  observer,  the  sun,  and  the  footprint: 

DThetaF  =  Pi/2.  -  Side (Pi/2. -DThetaO,  -DPsiPO,  DBeta) 

IF  (DPsiPO  .GE.  Pi)  THEN 

DPhiF  =  DPhiO  +  Angle(Pi/2. -DThetaO, DBeta, Pi/2. -DThetaF) 

ELSE 

DPhiF  =  DPhiO  -  Angle (Pi/ 2. -DThetaO, DBeta, Pi/2. -DThetaF) 

END  IF 

DTo  =  Side (Pi/2. -DThetaS,  DPhiS-DPhiF,  Pi/ 2 .-DThetaF) 

To  =  REAL (DTo) 

IF  (DPhiS  .GE.  DPhiO)  THEN 

DPo  =  DPsiW  +  Angle (Pi/2. -DThetaF,  Pi/ 2. -DThetaS,  DTo) 

ELSE 

DPo  =  DPsiW  -  Angle (Pi/2. -DThetaF,  Pi/ 2 .-DThetaS,  DTo) 

END  IF 

Po  =  REAL (DPo) 

IF  (DPhiF  .GE.  DPhiO)  THEN 

DPr  =  DPsiW  -  Angle (Pi/2. -DThetaF,  Pi/2 . -DthetaO,  DBeta) 

ELSE 

DPr  =  DPsiW  +  Angle (Pi/ 2. -DThetaF,  Pi/ 2 . -DthetaO,  DBeta) 

END  IF 

Pr  =  REAL (DPr) 

Calculate  specular  slope  (merely  for  print-out,  not  used  for 
further  calculations) ; 

Find  the  rectangular  receiver  coordinates, 

Ar  =  sin(Tr) *cos(Pr) 

Br  =  sin(Tr) *sin(Pr) 

Cr  =  cos(Tr) 

Find  the  the  rectangular  source  coordinates  for  the  solar  center, 
Ao  =  sin{To) *cos(Po) 

Bo  =  sin(To) *sin(Po) 

Co  ==  cos  (To) 

Find  the  Cox-Munk  wind  dependent  slope  variances 
if  (W  -ge.  .01)  then 

use  the  Cox-Munk  variance  for  a  real  sea 
Vu  =  0.000  +  3.16E-3*W 
Vc  =  0.003  +  1.92E-3*W 

else 

use  a  delta  function  for  an  ideal  calm  sea 
Vu  =  .0001 
Vc  =  .0001 

end  if 

pO  =  l./(2.*pi*sqrt(Vu*Vc)) 

Find  the  slopes  (Sxo,Syo)  for  a  specular  reflection  from 
source  (To,Po)  to  receiver  (Tr,Pr) , 
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if  (abs(Ao  +  Ar)  .le.  delta)  then 
Sxo  =0. 

else  if  ((Co  +  Cr)  .le.  delta)  then 
Sxo  =  sign (infinity,  -(Ao  +  Ar) ) 

else 

Sxo  =  -  (Ao  +  Ar)/(Co  +  Cr) 

end  if 

if  (abs(Bo  +  Br)  .le.  delta)  then 
Syo  =  0. 

else  if  ((Co  +  Cr)  -le.  delta)  then 
Syo  =  sign (infinity,  -(Bo  +  Br) ) 

else 

Syo  =  -  (Bo  +  Br)/(Co  +  Cr) 

end  if 

C  Calculate  the  Cox-Munk  tilt  and  slope: 
arg  =  Sxo**2/Vu  +  Syo**2/Vc 
p  =  p0*exp(-0.5*arg) 

Tn  =  atan(sqrt(Sxo**2  +  Syo**2)) 

C  and  print  to  "OUT": 

WRITE  (1P6,1000) 

WRITE  (IP6,1010)  DBeta*R2D,DPsiPO*R2D,AMOD(DPsiW*R2D, 360. ) 

IF  ((lEMSCTold)  .EQ.  2)  THEN 

WRITE  (IP6,1020)  DThetaO*R2D,DPhiO*R2D,DThetaF*R2D,DPhiF*R2D, 
+  DThetaS*R2D,DPhiS*R2D 

END  IF 

WRITE  (IP6,1030) 

WRITE  (IP6,1040)  Tr*R2D,AMOD(Pr*R2D,360. ) 

IF  ( (lEMSCTold)  .EQ.  2)  THEN 

WRITE  (IP6,1050)  To*R2D,AM0D(Po*R2D,360. ) ,Sr2d*Tn,sqrt(arg) ,p 
END  IF 
C 

1000  f ormat (/,' SUMMARY  OF  OBSERVATION  GEOMETRY'/) 

1010  format  (lOX, 'BETA  =',F10.5, 

+  '  DEG ' , / , 

+  lOX, 'PATH  AZIMUTH  =',F10.3, 

+  '  DEG  EAST  OF  NORTH',/, 

+  lOX, 'WIND  AZIMUTH  =',F10.3, 

+  '  DEG  EAST  OF  NORTH', \) 

1020  format  (1 OX, 'RECEIVER  LATITUDE  =',F10.3, 

+  '  NORTH  OF  EQUATOR',/, 

+  1 OX, 'RECEIVER  LONGITOTE  =',F10.3, 

+  '  WEST  OF  GREENWICH',/, 

+  1 OX, 'FOOTPRINT  LATITUDE  =',F10.3, 

+  '  NORTH  OF  EQUATOR',/, 

+  lOX, ' FOOTPRINT  LONGITUDE  = ' , FIO . 3 , 

+  '  WEST  OF  GREENWICH',/, 

+  lOX, 'SUBSOLAR  LATITUDE  =',F10.3, 

+  '  DEG  NORTH  OF  EQUATOR',/, 

+  lOX, 'SUBSOLAR  LONGITUDE  =',F10.3, 

+  '  DEG  WEST  OF  GREENWICH',//) 

1030  f ormat (/, 'VALUES  SEEN  FROM  FOOTPRINT'/) 

1040  format  (lOX, 'RECEIVER  ZENITH  ANGLE  =',F10.3, 

+  '  DEG',/, 

+  lOX, 'RECEIVER  AZIMUTH  =',F10.3, 


(A12) 
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+  *  DEG  WEST  OF  UP  WIND') 

1050  format  (lOX, 'SOLAR  ZENITH  ANGLE  =*,F10.3, 


-f 

'  DEG',/, 

lOX, 

' SOLAR  AZIMUTH 

=  1 

'  ,F10.3, 

'  DEG  WEST  OF  UP 

WIND* ,/, 

+ 

lOX, 

•SOLAR  SPECULAR 

TILT  = ’ 

\F10.3, 

'  DEG  ( ' ,  F6 . 2 , 

*  SIGMA, 

PROS  = ' , IPEIO 

*) 


return 

END 


SUBROUTINE  SunFoot (PsiO ,  DelO,  PsiPO,  Beta,  Psi) 

CU  USES  Angle,  Side 

COMMON  /Constants/  Spi, Sr 2d, Sd2r, epsilon, delta, onem, onep, infinity 
COMMON  /Geometry/  To,Po,Tr,Pr 

COMMON  /Card2/  IHAZE,  ISEASN,  IVULCN,  ICSTL,  ICLD,  IVSA,  VIS,W,  WHH, 

+  RAINRT 

COMMON  /IFIL/  IRD,IPR,IPU,NPR,IPR1,IP6,IP7,IP8,IP4,IRDS,IP6S, 

•f  ITR,Ipath,Is)cy,Isun 

COMMON  /Sea/  Sea,  Hit,  Msea,  TBOUNDold,  lEMSCTold 
REAL  infinity 


*********************************************************************** 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


This  routine  calculates  zenith  angles  and  azimuthal  angles  from 
the  footprint,  defined  as  the  location  where  a  hit  has  just 
occurred,  whenever  the  sun  is  involved.  (PsiPO  is  not  used  in 
the  calculation;  it's  passed  in  only  to  be  printed  out.) 

All  the  arguments  are  inputs,  and  all  are  angles  in  degrees: 

PsiO  is  the  solar  azimuth  measured  from  the  observer's 
line-of-sight  (+E  of  N) . 

DelO  is  the  solar  zenith  angle  as  seen  by  the  observer . 

Beta  is  the  angle  subtended  at  the  center  of  the  earth 
between  the  observer  and  the  footprint. 

Psi  is  the  wind  azimuth  measured  from  the  observer's 
line-of-sight  {-fE  of  N)  . 

(Psi  is  ASSUMED  to  be  the  same  at  the  footprint.) 

The  outputs  are  To,  Po,  and  Pr,  angles  in  radians  passed 
through  the  common  block  "Geometry" : 

To  is  the  zenith  angle  of  the  center  of  the  sun 
from  the  footprint- 

Po  is  the  azimuth  (+  W  of  PsiW)  of  the  center  of  the  sun 
from  the  footprint. 

Pr  is  the  azimuth  (+  W  of  PsiW)  receiver  as  seen 
from  the  footprint. 

Note:  Tr  has  been  calculated  in  DPFNMN  and  is  used 

in  this  subroutine  only  for  printing  to  "OUT”. 

Last  revision:  June  14,  1995. 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 
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************************************************************************ 


DOUBLE  PRECISION  DPsiO,  DDelO,  DBeta,  Pi,  Side,  Angle,  DTo,  DPo 
Pi  =  4.*DATAN{1.) 
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D2R  =  Pi/180. 

R2D  =  REAL (180. /Pi) 

First,  convert  to  radians  and  increase  precision: 

DPsiO  =  DBLE{PsiO)*D2R 
DDelO  =  DBLE(Del0)*D2R 
DBeta  =  DBLE(Beta)*D2R 
DPsi  =  DBLE(Psi)*D2R 

then  use  the  geometry  of  the  spherical  triangle  connecting 
the  observer,  the  sun,  and  the  footprint: 

DTo  =  Side (DBeta,  DPsiO,  DDelO) 

To  =  REAL (DTo) 

DPr  =‘Pi  +  DPsi 
Pr  =  REAL (DPr) 

If  (DPsiO  .GT.  0.)  then 

DPo  =  DPr  +  Angle (DTo,  DDelO,  DBeta) 

Else  if  (DPsiO  .EQ.  0.)  then 
DPo  =  DPr  +  Pi 

Else 

DPo  =  DPr  -  Angle (DTo,  DDelO,  DBeta) 

End  If 

Po  =  REAL (DPo) 

Calculate  specular  slope  (calculations  from  now  on  merely  for 
print-out,  not  for  further  calculations) : 

Find  the  rectangular  receiver  coordinates, 

Ar  =  sin(Tr) *cos (Pr) 

Br  =  sin(Tr) *sin(Pr) 

Cr  =  cos(Tr) 

Find  the  the  rectangular  source  coordinates  for  the  solar  center, 
Ao  =  sin(To) *cos (Po) 

Bo  =  sin(To) *sin(Po) 

Co  =  cos (To) 

Find  the  Cox-Munk  wind  dependent  slope  variances 
if  (W  .ge.  .01)  then 

use  the  Cox-Munk  variance  for  a  real  sea 
VU  =  0.000  +  3.16E-3*W 
Vc  =  0.003  +  1.92E-3*W 

else 

use  a  delta  function  for  an  ideal  calm  sea 
Vu  =  .0001 
Vc  ==  .0001 

end  if 

pO  =  1. / (2 • *pi*sqrt (Vu*Vc) ) 

Find  the  slopes  (Sxo,Syo)  for  a  specular  reflection  from 
source  (To,Po)  to  receiver  (Tr,Pr) , 
if  (abs(Ao  +  Ar)  .le.  delta)  then 
Sxo  =0. 
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else  if  ((Co  +  Cr)  .le.  delta)  then 
Sxo  =  sign  (infinity,  -(Ao  +  Ar) ) 

else 

Sxo  =  -  (Ao  +  Ar)/(Co  +  Cr) 

end  if 

if  (abs(Bo  +  Br)  .le.  delta)  then 
Syo  =  0. 

else  if  ((Co  +  Cr)  .le.  delta)  then 
Syo  =  sign  (infinity,  -(Bo  +  Br) ) 

else 

Syo  =  -  (Bo  +  Br)/(Co  +  Cr) 

end  if 

C  Calculate  the  Cox-Munk  tilt  and  slope: 

arg  =  Sxo**2/Vu  +  Syo**2/Vc 
p  =  p0*exp(-0. 5*arg) 

Tn  =  atan(sqrt (Sxo**2  +  Syo**2) ) 


C  and  print  to  *'0UT”: 


WRITE  (IP6,2000) 

WRITE  (IP6,2010)  DBeta*R2D,PsiP0,AM0D((Psi+PsiP0) ,  360.) 

WRITE  (IP6,2030) 

WRITE  (IP6,2  040)  Tr*R2D , AMOD (Pr*R2D, 360 . ) 

IF  ((lEMSCTold)  .EQ.  2)  THEN 

WRITE  (IP6,2050)  To*R2D, AMOD (Po*R2D, 360 . ) , Sr2d*Tn, sqrt (arg) , p 
END  IF 


^2000  fonnat(/,  'SUMMARY  OF  OBSERVATION  GEOMETRY*/) 


2010  format  (lOX, 
4- 

+  lOX, 

.+ 

+  lOX, 


BETA  =',F10.5, 

DEG ' , / , 

PATH  AZIMUTH  =’,F10.3, 

DEG  EAST  OF  NORTH',/, 

WIND  AZIMUTH  =',F10.3, 

DEG  EAST  OF  NORTH ' , \ ) 

2030  formate/, 'VALUES  SEEN  FROM  FOOTPRINT'/) 

■RECEIVER  ZENITH  ANGLE  =',F10.3, 

DEG',/, 

RECEIVER  AZIMUTH  = ' , FIO . 3 , 

DEG  WEST  OF  UP  WIND ' ) 

SOLAR  ZENITH  ANGLE  =',F10.3, 

DEG',/, 

SOLAR  AZIMUTH  = ' , FIO . 3 , 

DEG  WEST  OF  UP  WIND ' , / , 

SOLAR  SPECULAR  TILT  =',F10.3, 

DEG  (',  F6.2,  '  SIGMA,  PROB  =' , IPEIO. 3 , ' )  ' ) 


i  format 

(lOX, 

+ 

lOX, 

+ 

1  format 

(lox. 

+ 

lOX, 

+ 

lOX, 

+ 

return 

END 

SUBROUTINE  Card 

COMMON  /Card2/  IHAZE, ISEASN, IVULCN, ICSTL, ICLD, IVSA, VIS, WSS, WHH, 
+  RAINRT 

COMMON  /Card3/  HI, H2, ANGLE, RANGE, BETA, RE, LEN,Psi,SeaSwitch 
COMMON  /Card3Al/  IPARM, IPH, IDAY, ISOURC 

COMMON  /Card3A2/  PARM1,PARM2,PARM3,PAEM4,GMT,PSIP0,ANGLEM,G 
COMMON  /IFIL/  IRD,IPR,IPU,NPR,IPR1,IP6,IP7,IP8,IP4,IRDS,IP6S, 
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+  ITR, Ipath, Isky , I sun 

COMMON  /Constants/  pi, r2d,d2r, epsilon, delta, onem^onep,  inf inity 
COMMON  /Geometry/  To,Po,Tr,Pr 

COMMON  /Sea/  Sea,Hit  ,Msea,TBOUNDold,  lEMSCTold 
REAL  infinity 
LOGICAL  SeaSwitch 


************************************************************************ 


* 

*  Issues  new  MODTRAN  cards  for  the  sea  routines. 

* 

*  When  lEMSCTold  =  1,  no  sun  is  involved,  and  three  new  sky 

*  cards  are  issued  to  "TAPES. SEA":  one  for  Tmin,  the  minimum  sky 

*  zenith  angle  expected  at  the  current  wind  speed,  one  for 

*  for  Tmax,  the  maximum  zenith  angle  expected,  and  one  for  Tav, 

*  the  sky  zenith  angle  halfway  between  Tmax  and  Tmin. 

* 

*  When  lEMSCTold  =  2,  the  sun  is  involved,  and  after  each  new  sky 

*  card  the  original  cards  3A1  and  3A2  are  reissued.  At  the  very 

*  end  of  the  file  there  is  one  sun  card.  Hence  the  number  of  new 

*  cards  issued  to  'TAPES, SEA"  is  10  when  lEMSCTold  =  2. 

* 

*  Last  revised:  February  28,  1995. 

* 

*********************************************************************** 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


Irpt  =  3 

C  First,  find  the  wind~dependent  sky  angles  Tmin  and  Tmax: 

if  (WSS  .ge,  ,01)  then 

C  use  the  Cox-Munk  standard  deviation  for  a  real  sea 

Su  =  sqrt(3.16E-3*WSS) 

Sc  =  sqrt(3.E-3  +  1.92E-3*WSS) 

else 

C  use  a  delta  function  for  an  ideal  calm  sea 

Su  =  .01 
Sc  =  .01 

end  if 
S  =  2.8 

is  the  number  of  standard  deviations  to  which  the 
wave  slope  integral  will  be  carried;  for  S  =  2.8 
99  %  of  the  volume  under  the  distribution  is  captured. 
dT  =  2 , 02* (atan{S*amaxl (Su, Sc) ) ) 

Tmin  =  amaxl(Tr  -  dT,  1.) 

Tmax  =  aminl(Tr  +  dT,  d2r*89.) 

C  Next,  open  TAPES. SEA,  the  alternate  file  to  TAPES: 

open  (Irds,  file  =  'Tapes. Sea* ,  status  =  'unknown') 

C  then  write  the  sky  cards  (lEMSCT  =  2,  ITYPE  ==  3)  : 
do  Ts  =  Tmin,  Tmax,  onem* (Tmax-Tmin) /2 . 
write  (Irds,  150)  Irpt 

write  (Irds,  100)  0. ,0. ,Ts*r2d,0. ,0. ,0. , 0, Psi, SeaSwitch 
if  (lEMSCTold  .eq,  2)  then 

write  (Irds,  400)  IPARM,  IPH,  IDAY,  ISOURC 
write  (Irds,  500)  PARM1,PARM2  ,PARM3  ,PARM4  ,GMT,PSIPO, 
+  ANGLEM,G 
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end  if 
end  do 


C  write  the  sun  card  (lEMSCT  =  3,  ITYPE  =  3)  if  necessary: 
if  (lEMSCTold  .EQ.  2)  then 
write  (Irds,  150)  Irpt 

write  (Irds,  200)  0. , 0. ,To*r2d, IDAY, 0. , 0, 0. 
end  if 

C  and  finally 
rewind  Irds 

C  so  it  can  be  read  from  the  beginning  by  the  driver. 


return 


100  format 
150  format 
200  format 
400  format 
500  format 


(6F10.3,I5,F10.3,L5) 

(15) 

(3F10.3,I5,5X,F10.3,I5,F10,3) 

(415) 

(8F10.3) 


END 


***********  ********  FUNCTIONS  ******************************** 

*  * 

*  Latest  revision;  May  5,  1994  for  Side  and  Angle.  * 

*  .  * 
************************************************************************ 

FUNCTION  Side(a,C,b) 

is  the  Law  of  Cosines  for  a  spherical  triangle  with  sides  a,  b, 
and  c  and  opposite  angles  A,  B,  C.  The  three  parameters  are 
the  two  sides  a  and  b  and  the  included  angle  C.  Side  is  the 
value  of  side  c  opposite  the  included  angle  C.  Angles  are  in 
radians. 

double  precision  a,C,b,Side 

Side  =  dacos(dcos (a) *dcos (b)  +  dsin (a) *dsin(b) *dcos (C) ) 

END 

FUNCTION  Angle(a,c,b) 

is  the  Law  of  Cosines  for  a  spherical  triangle  with  sides  a,  b, 
and  c  and  opposite  angles  A,  B,  C.  The  three  parameters  are 
the  three  sides.  Angle  is  the  value  of  angle  C  opposite  side 
c,  the  middle  parameter  in  the  list.  Angles  are  in  radians, 

double  precision  pi,a,b,c,Angle,Arg,aa,bb 
pi  =  4.*datan(l.) 
aa  =  dminl(a,  b) 
bb  =  dmaxl(a,  b) 

if  (abs(aa)  .le.  l.D-5)  then 

Angle  =  dacos(dtan(aa) /dtan(bb) ) 
if  (abs(bb)  .It.  abs(c))  Angle  =  pi  -  Angle 

end  if 
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Arg  =  (dcos(c)  -  dcos(a) *dcos(b) ) / (dsin(a)*dsin(b) ) 
if  (abs(Arg)  .ge.  (1  -  l.D-14))  then 
Angle  =  0. 

else 

Angle  =  dacos(Arg) 

end  if 


FUNCTION  Rho{ Omega,  V) 

CU  USES  SeaData 

COMMON  /Sealndex/  AlphaOl (100)  ,  Alpha02(20), 

+  BetaOl  (100),  Beta02  (20) 

************************************************************************ 

* 

* 

*  Calculates  reflectivity  of  sea  water  between  52.63  cm-l  and  * 

*  25,000  cm-l  using  equations  (74)  and  (78)  of  Stratton,  "Electro-  * 

*  magnetic  Theory",  1941,  page  505,  ff.  The  sea  water  is  assumed  to  * 

*  be  a  conducting  medium;  both  real  and  imaginary  parts  of  the  * 

*  index  of  water  are  used.  The  notation  of  Stratton  is  adhered  to  * 

*  as  closely  as  possible,  ^ 

* 

*  Omega  is  the  angle  of  incidence  in  radians;  V  is  the  wave-  * 

*  number  in  cm-1;  Rho  is  the  reflectivity.  * 

* 

*  Last  revision:  November  28,  1994  * 

* 

* 

************************************************************************ 

C  The  four-point  interpolation  functions  are: 

WM(X)  -  (X  -  l.)*(X  -  2.)*X/6. 

W0(X)  -  (X  -  l.)*(X  -  2.)*(X  +  l.)/2. 

W1(X)  =  (X  +  l.)*(X  -  2.)*X/2. 

W2(X)  =  (X  +  10*{X  -  l,)*X/6. 

IF  ((Omega  .LT.  0.)  .OR.  (Omega  .GT.  1.57080))  THEN 
C  Omega  is  out  of  bounds;  set  reflectivity  to  0%  and  return: 

Rho  =  0. 

RETURN 
END  IF 

IF  (V  .EQ.  0.)  THEN 

C  set  reflectivity  to  100%  and  return: 

Rho  =  1. 

RETURN 
END  IF 


W  =  1.E4/V 


IF  (W  .LT.  0.399999)  THEN 
C  print  error  message  and  return: 

Rho  =0. 

WRITE  (6,  1000)  V 
RETURN 
END  IF 
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IF  (0.4  .LE.  W  .AND.  W  .LE.  19.8)  THEN 
C  use  Lagrange  4  point  interpolation  on  Block  01  data  which 

C  are  at  0.2  um  spacing  between  0.2  and  20  um: 

I  =  INT(W/0.2) 

Fr  =  M0D(W,  0.2) /O. 2 

Alphal  =  W2(Fr)*Alpha01(I  +  2)  -  W1 (Fr) *Alpha01 (I  +  1) 

+  +  WO(Fr)*Alpha01(I)  -  WM(Fr) *Alpha01(I  -  1) 

Betal  =  W2(Fr)*Beta01  (I  +  2)  -  W1 (Fr) *Beta01  (I  +  1) 

+  +  W0(Fr)*Beta01  (I)  -  WM(Fr) *Beta01  (I  -  1) 

END  IF 

IF  (19.8  .LT.  W  .AND.  W  .LT.  190.)  THEN 
C  use  Lagrange  4  point  interpolation  on  Block  02  data  which 

C  are  at  10  um  spacing  between  10  and  200  um: 

I  =  INT(W/10.) 

Fr  =  MOD(W,  10. ) /lO. 

Alphal  =  W2(Fr)*Alpha02(I  +  2)  -  W1 (Fr) *Alpha02 (I  +  1) 

+  +  W0(Fr)*Alpha02(I)  -  WM(Fr)*Alpha02(I  -  1) 

Betal  =  W2 (Fr) *Beta02  (1+2)  -  Wl(Fr) *Beta02  (I  +  1) 

+  +  WO(Fr) *Beta02  (I)  -  WM(Fr) *Beta02  (I  -  1) 

END  IF 

IF  (190.  .LE.  W)  THEN 
C  print  error  message  and  return: 

RHO  =  0. 

WRITE  (6,  1000)  V 
RETURN 
END  IF 

G  =  Alphal**2  -  Betal**2  -  SIN(Omega) **2 
H  =  4*(Alphal**2)*(Betal**2) 

P  =  SQRT(0.5*(-G  +  SQRT(H  +  G**2) ) ) 

Q  =  SQRT(0.5*(+G  +  SQRT(H  +  G**2))) 

C  Stratton,  Equation  (74) ,  p.  505: 

C  =  (Q  -  COS (Omega) )**2  +  P**2 

D  =  (Q  +  COS (Omega) ) **2  +  P**2 

Rp  =  C/D 

C  Stratton,  Equation  (77) ,  p.  506: 

E  =  ((Alphal**2  -  Betal**2)*COS (Omega)  -  Q)**2 

F  =  ((Alphal**2  -  Betal**2) *COS (Omega)  +  Q)**2 

T  =  (2*Alphal*Betal*COS (Omega)  -  P)**2 

U  =  (2*Alphal*Betal*C0S (Omega)  +  P) **2 

Rs  =  (E  +  T)/(F  +  U) 

Rho  =  (Rp  +  Rs)/2. 

RETURN 

1000  FORMAT  (•  *****  WARNING  -  INPUT  FREQUENCY  =  ',  1PG12.5,  'CM-l', 

+  /,  '  OUTSIDE  VALID  RANGE  OF  52.63  TO  25,000  CM-l  ******',/) 

END 

BLOCK  DATA  SeaData 

COMMON  /sealndex/  AlphaOl(lOO) ,  Alpha02(20), 

+  BetaOl  (100),  Beta02  (20) 
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************************************************************************ 
*  * 

*  These  data  for  the  optical  index  of  water  have  been  taken  from  * 

*  the  literature.  From  0.2  to  2  microns  (blocks  01  up  to  second  * 

*  entry  of  row  B)  and  from  10  to  200  microns  (blocks  02)  the  data  * 

*  are  from  G.  M.  Hale  and  M.  R.  Querry,  "Optical  Constants  of  Water* 

*  in  the  200-nm  to  200-um  Wavelength  Region,"  Appl.  Opt.  3,  555  * 

*  (1973).  These  data  are  for  pure  water.  * 

*  * 

*  From  2.2  to  20  microns  (blocks  01  from  the  third  entry  of  row  B  * 

*  to  the  end)  the  data  are  from  M.  R.  Querry,  W.  E.  Holland,  R.  C.  * 

*  Waring,  L.  M.  Earls,  and  M.  D.  Querry,  "Relative  Reflectance  and  * 

*  Complex  Refractive  Index  in  the  Infrared  for  Saline  Environmental* 

*  Waters,"  J.  Geophys.  Res.  82,  1425  (1977),  Table  3,  Pacific  * 

*  Ocean  columns.  These  data  are  for  salt  water,  * 

*  * 
************************************************************************ 


c 

c 


c 

c 


c 

c 


Real 

part  of 

the  index  of 

sea  water  from 

0.2  to 

20  microns  in 

steps 

of  0.2 

microns 

DATA  AlphaOl 

/ 

A 

1.396, 

1.339, 

1.332, 

1.329, 

1.327, 

1.324, 

1.321, 

1.317, 

B 

1.312, 

1.306, 

1.303, 

1.287, 

1.251, 

1.151, 

1.384, 

1.479, 

C 

1.421, 

1.388, 

1.368, 

1.355, 

1.347, 

1.339, 

1.335, 

1.335, 

D 

1.332, 

1.324, 

1.312, 

1.296, 

1.268, 

1.271, 

1.371, 

1.353, 

E 

1.340, 

1.330, 

1.324, 

1.319, 

1.314, 

1-307, 

1-302, 

1-297, 

F 

1.291, 

1.286, 

1.279, 

1.272, 

1.268, 

1.264, 

1.258, 

1.249, 

G 

1.240, 

1.229, 

1.218, 

1.204, 

1.190, 

1,177, 

1.165, 

1.151, 

H 

1.140, 

1.132, 

1.124, 

1.119, 

1.121, 

1.126, 

1.134, 

1.142, 

I 

1.152, 

1.164, 

1.177, 

1.189, 

1.201, 

1.212, 

1.224, 

1.234, 

J 

1.242, 

1.253, 

1.261, 

1.273, 

1.284, 

1.296, 

1.309, 

1.320, 

K 

1,331, 

1.339, 

1.349, 

1.358, 

1.366, 

1.379, 

1.390, 

1.399, 

L 

1.408, 

1.417, 

1.426, 

1.435, 

1.443, 

1.450, 

1.458, 

1.464, 

M 

1.470, 

1.474, 

1.477, 

1.480 

/ 

Real  part  of  the  index  of  sea  water  from  10  to  200  microns  in 
steps  of  10  microns: 

DATA  Alpha02  / 

A  1.218,  1.480,  1.551,  1.519,  1.587,  1.703,  1.821,  1.886, 

B  1.924,  1.957,  1.966,  2.004,  2.036,  2,056,  2.069,  2.081, 

C  2.094,  2.107,  2.119,  2.130  / 


Imaginary  part  of  the  index  of  sea  water 
steps  of  0.2  microns: 

DATA  Beta 01  / 


A 

0.000, 

0.000, 

0.000, 

0.000, 

0.000 

B 

0.000, 

0.001, 

0.000, 

0.001, 

0.003 

C 

0.018, 

0.005, 

0.003, 

0.004, 

0.007 

D 

0.013, 

0.011, 

0.010, 

0.013, 

0.032 

E 

0.035, 

0.033, 

0.032, 

0.031, 

0,031 

F 

0.034, 

0.036, 

0.038, 

0.041, 

0.044 

G 

0.048, 

0.050, 

0.054, 

0.060, 

0.068 

H 

0.125, 

0.145, 

0.166, 

0.191, 

0.216 

I 

0.297, 

0.313, 

0.326, 

0.338, 

0.347 

J 

0.377, 

0,385, 

0.393, 

0.401, 

0.407 

K 

0,420, 

0.422, 

0.424, 

0.427, 

0.430 

L 

0.431, 

0.430, 

0.429, 

0.427, 

0,425 

from  0.2  to  20  microns  in 


0.000, 

0.000 

0.114, 

0.263 

0,011, 

0.016 

0.108, 

0.087 

0.032, 

0.032 

0.046, 

0.046 

0.079, 

0.091 

0.239, 

0.260 

0.357, 

0,363 

0.413, 

0.417 

0.432, 

0.432 

0.423, 

0.420 

0.000, 

0.085, 

0.016, 

0.044, 

0.033, 

0.047, 

0.107, 

0.279, 

0.371, 

0.418, 

0.432, 

0.416, 
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M  0.414,  0.410,  0.406,  0.393  / 

C  Imaginary  part  of  the  index  of  sea  water  from  10  to  200  microns  in 
C  steps  of  10  microns: 

DATA  Beta02  / 

0.051,  0.393,  0.328,  0.385,  0.514,  0.587,  0.576,  0.547, 
0.536,  0.532,  0.531,  0.526,  0.514,  0.500,  0.495,  0.496, 
0.497,  0.499,  0.501,  0.504  / 

END 
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Washington,  DC  20375  (3) 

Naval  Surface  Warfare  Center 
Bethesda,  MD  20084-5000 

Naval  Air  Warfare  Center, 

Weapons  Division 

Point  Mugu,  CA  93042-5001  (3) 


Naval  Postgraduate  School 

Monterey,  CA  93943-5 100  (10) 

ARETE'  Associates 

Sherman  Oaks,  CA  91413  (3) 

Photon  Research  Associates,  Inc. 

La  Jolla,  CA  92037-1020 

Viz  Ability,  Inc. 

Corvallis,  OR  97330 

Michigan  Technical  University 
Houghton,  MI  4993 1-1295  (4) 

University  of  California,  San  Diego 
La  Jolla,  CA  92093-0230  (2) 


Office  of  Naval  Research 

Arlington,  VA  22217-5660  (4) 

Aerodyne  Research,  Inc. 

Billerica,  MA  01821-3976  (2) 

Science  &  Technology  Corporation 
Hampton,  VA  23666-1340 

Jet  Propulsion  Laboratory 
Pasadena,  CA  91109-8099 

Texas  A&M  University 
College  Station,  TX  77843 
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FOREIGN 

Physics  &  Electronics  Laboratory  TNO 
The  Netherlands  (4) 

Physics  Department  UMIST 
United  Kingdom 

European  Centre  for  Medium  Range  Weather  Forcasts 
United  Kingdon 

Defense  Research  Establishment  Valcartier 
Quebec  GOA  IRO  Canada  (3) 


